FINAL REPORT 


ON 

NASA PROJECT NAG 3-1376 


DEVELOPMENT OF A MODEL BASED TECHNIQUE 

FOR 

GEAR DIAGNOSTICS 
USING THE WIGNER-VILLE METHOD 


by 

F. Choy 
A. Xu 

V. Polyshchuk 

Department of Mechanical Engineering 
The University of Akron 
Akron, OH 44325-3903 



ABSTRACT 


Imperfections in gear tooth geometiy often results from errors in the 
manufacturing process or excessive material wear during operation. Such faults in the 
gear tooth geometry can result in large vibrations in the transmission system, and, in some 
cases, may lead to early failure of the gear transmission system. This report presents the 
study of the effects of imperfection in gear tooth geometiy on the dynamic characteristics 
of a gear transmission system. The faults in the gear tooth geometry are modeled 
numerically as the deviation of the tooth profile from its original involute geometiy. The 
changes in gear mesh stiffness due to various profile and pattern variations are evaluated 
numerically. The resulting changes in the mesh stiffness are incorporated into a computer 
code to simulate the dynamics of the gear transmission system. A parametric study is 
performed to examine the sensitivity of gear tooth geometry imperfections on the 
vibration of a gear transmission system. The parameters varies in this study aconsist of 
the magnitude of the inperfection, the pattern of the profile variation, and the total 
number of teeth affected. Numerical results from the dynamic simulations are examined 
in both the time and the frequency domains. A joint time-frequency analysis procedure 
using the Wigner-Ville Distribution is also introduced to identify the location of the 
damaged tooth from the vibration signature. Numerical simulations of the system 
dyanmics with gear faults were compared to experimental results. An optimal tracker was 
introduced to quantify the level of damage in the gear mesh system. Conclusions are 
drawn from the results of this numerical study. 
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CHAPTER 1 


INTRODUCTION 


Gearing, one of the most universally used machine elements, is applied in 
mechanical systems of every size and description: fiom the tiny pinions in a watch or a 
computer system to the high-speed, heavily-loaded reduction gears of an aircraft gas 
turbine. In the last two decades, the use of gear transmissions in both defense and 
commercial applications has substantially increased. With the demand for higher power 
and performance, premature failures in transmissions often result in financial losses, and 
sometimes even lead to catastrophic consequences. In the aerospace industry, where both 
weight-to-load factor and efficiency are pushed to their design limits, one of the major 
concerns is fatigue failures in rotorcraft gear transmission systems. Such failures are 
often a result of excessive gear tooth wear and tooth crack formation. Presently, the 

prevention and management of premature equipment Mures has become a vital part of 
the maintenance program. 

Current on-board condition monitoring systems for gas-turbine engine systems 
often fail to provide sufficient time between warning and failure such that safety 
procedures can be implemented. On the other hand, inaccurate interpretation of 
operational conditions may result in false alarms and unnecessary repairs and downtime. 
The early detection of incipient failure in a mechanical system is of great practical 
importance as it permits scheduled inspections without costly shutdowns, along with 
indication of urgency and location for repair before any catastrophic Mure. 

Presently, a considerable amount of work in machine life prediction has been 
carried out using machine reliability and design life approaches[ 1,2 ]. However, most of 
this work is based on statistical predictions developed by Lundberg and Palmgren[ 3 ] 
without considering the conditions of machine components during various phases of their 
lifespan. Besides the work reported in recent years by the Mechanical Failure Prevention 
Group(MFPG)[ 4 - 6 ], very little has been cited in the literature concernin g condition- 
based Mure prediction. 
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The increasing requirements for long life and safe operation in mechanical 
systems call for the development of an accurate machine fault identification and failure 
prognostication system which is capable of on-line machine health monitoring as well as 
machine life prediction. One of the advanced fault identification procedures used in 
rotorcraft mechanical systems is the signature analysis of machine vibration/acoustic 
signals! 7 - 9 ]. The acquired machine vibration/acoustic signature is compared to a data 
bank of standard healthy machine operation signatures to pinpoint the abnormalities of 
the input signal. This procedure does not require a shut down of the rotating machinery, 
and can be used as an in-flight diagnostic and trend monitoring device. 

In order to develop an accurate machine fault identification and failure 
prognostication system, it is very important to understand the dynamics of a gear 
transmission system under a variety of fault conditions, as well as under nominal 
conditions. This study deals only with the basic knowledge in fault identification of gear 
component and the dynamic modelling of a gear transmission system under the effects of 
gear tooth geometry imperfection. 

The major objective of the research presented herein is the development of a 
numerical model to predict the vibration in a gear transmission system due to gear tooth 
imperfection or damage. The imperfection/damage of the gear tooth geometry is modeled 
numerically as the deviation of tooth profile from its original geometry. To simulate gear 
failures, a computer code developed at NASA Lewis Research Center! 10 ] was modified 
to simulate various types of gear mesh conditions. The changes in mesh stiffness due to 
the effects of tooth wear can be represented by using a tooth profile modification 
procedure. The resulting changes in mesh stiffness are incorporated into the dynamic 
simulation of the gear transmission system. 

In simulating the vibration of the transmission system, the equations of motion are 
established individually for each rotor-gear-bearing system. These localized changes in 
the gear mesh stiffness are incorporated into each gear-rotor model during the dynamic 
simulation [11-13]. The dynamics of each gear-rotor system are coupled with each other 
through the gear mesh interacting forces and the bearing supporting forces. The global 
vibrations of the system are evaluated by solving the set of coupled transient dynamics 
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equations of all the rotor systems simultaneously including the vibration of the casing. In 
order to minimize the computational effort, the number of degrees-of-freedom of the 
system are reduced by using a modal synthesis procedure[ 1 1,12 ]. 

Results from the model were examined using a joint time-frequency analysis 
method. This approach was chosen because the joint time-frequency analysis will 
provide an instantaneous frequency spectrum of the system at every instant of the 
revolution of the pinion while the traditional frequency analysis can only provide the 
average vibration spectrum of the signal. In other words, the time-changing spectral 
density from the joint time-frequency spectra will provide vital information concerning 
the frequency distribution concentrated at a particular instant. The Wigner-Ville 
Distribution (WVD) [ 14-16 ] was chosen for the joint-time frequency analysis. 
Considerable success has been achieved in applying the WVD to gear transmission 
systems [ 17,18 ] to recognize faults at various locations of die gear. 

In addition, experimental results obtained from a gear failure test rig at NASA 
Lewis Research Center were also used to experimentally validated the identification 
procedure using the WVD as well as to verify the numerical simulations. A technique for 
quantifying the damage in a gear mesh using an optimal tracker was also developed, and 
results obtained using the optimal tracker were compared with those from the 
experimental studies. 

Based on results of this study, general conclusions are made concerning the 
effects of local damage on the global damage of the gear transmission system. This study 
applied the above discussed methodology on a variety of damaged gear models. The 
numerical model used to simulate the dynamics of a gear transmission system with gear 
tooth geometry imperfections was successfully developed. In addition, considerable 
success was achieved in generating a comprehensive database of the vibration signal due 
to various kinds of gear tooth geometry imperfection patterns a in gear transmission 
system. Some limited success was achieved in quantifying fire damage using an optimal 
tracker. Using the developed analytical/numerical model, a gear more extensive 
fault/damage database can be developed for machine fault identification and failure 
prognostication research. 



CHAPTER 2 


THE EFFECTS OF GEAR TOOTH IMPERFECTION 
ON GEAR MESH STIFFNESS 


2.1, Objective 

The magnitude variation in the tooth mesh stiffness can affect gear mesh 
dynamics and loading significantly. The first step in determining the effects of gear tooth 
imperfections on mech dynamics is to determine the relationship between gear tooth 
imperfections and the resulting change to gear mesh stiffness. 

The objective of this chapter is to develop a relationship describing the effects of 
gear tooth imperfections on the static behavior of a gear system, with an emphasis on the 
gear mesh stiffness. The method which calculates the tooth mesh stiffness of a gear 
system is presented by Richardson [ 19 ]. The imperfection of a gear tooth is modeled 
numerically as the deviation of the tooth profile from its designed geometry. The changes 
in gear mesh stiffness due to various profile changes and imperfection patterns are 
evaluated numerically. A computer code developed at NASA Lewis Research Center[ 10 
] was modified to simulate various types of gear mesh conditions. 

2.2. Gear Kinematic Properties 

The ideal kinematic requirement for gear action is constant speed ratio. That is, 
the angular velocity of the driven gear should be a constant multiple of the angular 
velocity of the driving gear. Two curves that possess fee property of constant speed ratio 
when operated as contacting tooth surfaces are called conjugate curves. From fee infini te 



m 
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variety of possible conjugate curves, the involute has been almost universally accepted 
for use in gearing. An involute curve is generated by the end of a line that is unwound 
from the circumference of a circle called the base circle. 

Figure 1 shows two gear teeth in contact. Point L on gear 2 is in contact with 
point M on gear 1 . At this point of contact, the two tooth surfaces must be tangent to 
each other and consequently must have a common normal W 1 W 2 passing through the 
point of contact which intersects the line of centers Q C 2 at the instant center P, called the 
pitch point. Since ideal gears are assumed rigid, the velocities of L and M along the 
normal WiW 2 must be equal. The velocities ofL and M perpendicular to the normal 
sliding or relative velocity of the tooth surfaces. 

V l2 =w 2 C 2 W 2 ; J (2.1) 

where wj and w 2 are the angular velocities of gears 1 and 2, respectively. Hence, by 
similar triangles 


w, R, 


( 2 . 2 ) 


This equation is frequently used to define the law of gearing, which states that the pitch 
point must remain fixed on the line of centers. This means that all the lines of action for 
every instantaneous point of contact must pass through the pitch point. Consequently, 
ideal gears can be represented kinematically by two imaginary cylinders of radii Ri and 
R 2 , called pitch cylinders, which roll on each other without slipping . 




Figure 1 . Involute gear geometry 
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If no friction is present between the mating gear profiles, then the resultant force 
transmitted at the contact point L must lie along the common normal. For this reason the 
common normal is called the pressure line , and the angle between the normal and a line 
perpendicular to the line of centers C 1 C 2 is called the pressure angle 9. The locus formed 
by all points of contact as the gears rotate is known as the path of contact. 

In order to maintain continuous conjugate action, a series of conjugate curves are 
spaced uniformly around the circumference of a gear. The separation of these curves, 
measured along the pitch circle, is called the circular pitch 



(2.3) 


where i is the number of teeth, and R is the pitch-circle radius. The diametrical pitch is 
defined as the number of teeth on the gear per inch of pitch diameter, as indicated in the 
following equation; 


p -Tr < 24 > 

The relationship between circular pitch and diametral pitch is as follows: 

P c P = k (2.5) 

In an involute gear, the spacing of successive involutes along the pressure line or 
line of action is known as the normal pitch, and is related to the circular pitch defined by 
Equation (2.3) in the following way 

P n =P c cos6 (2.6) 


The contact-ratio is defined as the path of contact divided by the normal pitch, 
and is a measure of the average number of tooth-pairs in contact To provide continuous 
action the contact ratio must be greater that one, and for most power transmi ss ion 



gearing, the value of this quantity lies between the values of one and two. In some 
instances, the contact ratio can be as high as four. 

The radial length of the teeth beyond the pitch circle is called the addendum 
distance, and the radial depth of the teeth below the pitch circle is called the dedendum 
distance. By trade association standards, these distances are specified as constant 
multiples of the circular pitch. 

The location in the gear mesh of the contact point between mating teeth can be 
specified conveniently by the distance, s, between the pitch point and the contact point, 
measured along the line of action. 

When load is being transmitted through the gear mesh, the load is carried either by 
one pair of teeth alone or jointly by two or more pairs of teeth depending on the value of 
the contact ratio. It is assumed here that the contact ratio is between one and two. As the 
gears rotate, the load is transferred from teeth that are in mesh to succeeding teeth that are 
moving into mesh. Similarly, teeth moving out of mesh relinquish load as they leave 
contact. 

2.3. Load-Deflection Properties of a Gear Mesh 
2.3. 1 Definition of Spring-Stiffness of Gear Mesh 

The ideal curves that are used to form gear teeth are designed to produce a 
constant speed ratio. That is , the gears behave like two im aginar y pitch cylinders which 
roll without slipping. 



9 


In actual gears, the materials employed cannot be absolutely rigid; consequently, 
the gear-teeth will deflect due to the transmitted loads, which will cause the ideal pitch 
circles to slip. Thus a deviation from ideal kinematic operation occurs. 

For example, in Figure 2, gear 2 is fixed, by definition its pitch circle is also 
fixed. Now consider a torque ti to be applied to the mating gear 1. This torque on gear 1 
must be balanced, for static operation, by the moment of the resultant force F, which, in 
the absence of friction, acts along the pressure line. 

r,=Fcos0R, (2.7) 

or in terms of T, the component of W which acts tangentially to the pitch circle, 

ti = TR! (2.8) 

When friction is present, or contact between mating teeth lies off the pressure line, W and 
T in Equations. (2.7) and (2.8) no longer represent tooth loads exactly, but are still 
convenient ways of expressing the input torque tj. 

The spring stiffness of the gear mesh is defined as the amount of tangential load 
T, computed from Equation (2.8), to produce one unit of pitch-circle slip, 8, as shown in 
Figure 2 



(2.9) 
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Figure 2 (a) Pitch - Circle Slip 
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Figure 2 Definition of Gear-Mesh Spring-Stiffness 



12 


This definition can be equally well stated as the amount of load W acting along 
the pressure line, required to produce one unit of relative displacement(sr) between gears, 
measured along the pressure line 

k = ^ (2.10) 

Where 

F cos 0 = T ( 2 .H) 

and 


Sf = 5 cos 0 

These two spring stiffness are related by virtue of Equations. (2.1 1) and (2.12) 

T K 0 


( 2 . 12 ) 


k = 


S( cos 2 #) cos 2 6 


(2.13) 


2.3.2 Compliance for a Single Pair of Teeth 

The determination of the compliance of gear teeth is considerably difficult 
because it is an integral function of the entire loaded tooth. In addition, because of the 
stubbiness of the teeth, the foundation and shear effects are important. Cornell’s method 
[ 20 ] parallels, to a great extent, Weber’s work[ 21], however Cornell uses O’Donnell’s 
foundation factors [22,23 ]. The total compliance or flexibility of a gear tooth at the point 
of load, yr, is made up of three deflections: 

1 . The basic deflection of the tooth as a cantilever beam, y^; 

2. The deflection of the tooth caused by the fillet and foundation flexibility [ 22 ], yp; 

3. The local deflection caused by the contact and compression between the two teeth, yt. 
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When two gear teeth are in contact, the total compliance is their combined 
deflection per unit of load at the contact position or 

c =(yn + y T2 )/L = [(y B{ + y B2 ) + {y FX + y F2 ) + y L )]/L (2.14) 

where L is the applied load, and the deflections y are in the direction of the load. The 
methods to calculate the three types of deflections in equation (2.14) are given below. 


2.3.2. 1 The Basic Deflection of the Tooth as a Cantilever Beam 

When a load L is applied to the surface of a gear tooth, a deflection of the tooth 
occurs in the direction of the load. Suppose that the tooth is rigid near the point of 
loading. Then deflections of the tooth will still occur due to each of the following effects. 

1 . Bending of the tooth in the manner of a cantilever beam. 

2. Direct compression of the tooth due to the radial component of the load (N) 

3. Direct shearing of the tooth due to the tangential component of the load (T). 

The deflection and , therefore, compliance of a gear tooth over its beam portion is 
easily obtained using elementary strength of materials. Referring to Figure 3(a), the total 
of the bending and shear deflection in the direction of and at the point of application of 
the tooth load L, which is at radius Rl or position S along the line of action, be 

expressed in integral form as 


Lcos 2 <f>' L f I, rj 2.40(1 + //) + Tan 2 <}>[ , 

y * - e 1 i d7]+ 1 }dz 


(2.15) 


where I is the section modulus of the tooth, and using the relationship of 

5 = x or q=z; jj = (/"-£) =z = (/-*) (2.16) 
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and a value of 1 .2 has been assumed for the shear factor based on a rectangular tooth; 

G=E/ 2(1 + //) (2.17) 

and A is the cross-sectional area as a function of x or z. The deflection of the basic tooth 
can also be defined in a summation expression rather than integral form, see Figure 3(b), 
In this case the tooth beam deflection at and in the direction of the load is 
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y B = 


L cos 2 <f>' L 


Is, 


I, -'A + 3 s! 


(2A(\ + n) + Tan 2 <f>' L ) 


t=\ 


A 


(2. 18) 


Where 

l/7,-(V/,+W.,)/2 (2.19) 

and 

lA=0M+lM +l )/2 (2.20) 

Using these inverse forms for the values of I t and A t improves the accuracy for a small 
number of elements. In equation (2.18) 

/,=(/-*,) ( 2 ^ 1 ) 

and 

A = ~*i) ( 2 . 22 ) 

Both approaches for beam flexibility assume a narrow tooth width, W. For wide teeth 
where W/hp > 5, the flexibility is decreased by the antielastic effect, so that the values of I 
in equations (2.15) and (2.18) should be divided by (1-p 2 ). 
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23.2.2 The Deflection of the Tooth Caused by the Fillet and Foundation Flexibility 
Because of the fillet and the flexibility of the material to which the tooth is 
attached, additional deflection will occur at the load; see references [21-23 ]. This fillet 
and foundation deflection in the direction of the load, yp, is a function of the fillet 
geometry and the load position and direction, and is determined by the effective fillet 
length or angle yp for which the maximum deflection or work occurs at the load. 

Based on Figure 4, O’Donnell [ 22, 23 ] shows that the deflection at and in the 
direction of the tooth load due to the foundation effects, ypp, for plane stress is 


_ LCqsV l 16.67 
WE n 



+2(i -46 


+ 1334 


1 + 


Tan 2 ^, ^ 
2.4(1 + /J 


(2.23) 



The O’Donnell coefficients in equations (2.23) and (2.24) differ slightly from those given 
by Weber[ 21]. The first term in the brackets is the deflection at L due to the rotation 
caused by the moment at h F . The second term is the sum of the deflections at L due to the 
displacement at hp caused by the moment at h F and the rotation at h F caused by the shear 
force at h F . The first part of the third term is the displacement at L due to the shear force 
at h F based on the assumption that the effective depth for dete rmining this deflection is 
2 !i times the tooth thickness [ 21 ]. The second part of the third term is the deflection at 
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L due to the normal component of the load assuming the same relationship holds as 
indicated by the beam deflection equation; see equations (2.15) and (2.18). 

Referring to Figure 4, the deflection at and in the direction of the load due to the 
flexibility of the fillet and foundation, yp, the shaded region, is obtained fiom the 
summation of the fillet beam deflection, y^, using equations (2.15) and (2.18), and the 
foundation deflection ypF, using equation (2.23) or (2.24), i.e., 

yF = yFF +yre (2.25) 

where 


l F = l +r(Siny F -Siny) 
h F = h + 2 r(Cosy - Cosy F ) 
h = l + r(Siny j - Siny) 

= h + 2r{Cosy — Cosy i ) 


(2.26) 


The value of yp is the one that maximizes the value ofyp or yr> which can easily be done 
as one integrates or sums up the deflection of the tooth starting at the beginning of the 
fillet or at the load position. 


2.3.2.3 The Local Deflection Caused by Compression of the Tooth Surfaces 

The local compliance, yt, consists of the Hertz or line contact deflection plus the 
compression of each tooth between the point of contact and the tooth centerline. Figure 
1 5 gives the nomenclature for the parameters that determine the local deformation. There 
are several approaches to calculate the local compliance. All of the expressions for the 
local deformation are nonlinear with load because of the Hertz half contact width b. 

Here, the closed form approach developed by Weber[ 21 ] was adopted. 
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Weber in referencef 21] developed an expression specifically for the local 
deformation of two gear teeth, based on Hertz’s work on deformation between cylinders. 
In order to obtain a closed form solution, he assumed small deformations so that just the 
first two terms of the binomial expansion of the deformation needed to be used - i.e.. 





Figure 4. Fillet and foundation compliance of a gear tooth 
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b = 


4L 


nW 


n\ 

l E 2 , 


'[l/ r, + l/r 2 ] 


1/2 


(2.28) 


If Ei - E 2 and pi - J 4 . 2 , equation (2.27) reduces to 


y L = 


4(1- // 2 ) z 


;r£ IF 


r /==->( 


Inl 


4*a 


.2(1 -mV 


(2.29) 


The overall compliance, C, of the tooth pair is obtained by adding the local gear 
tooth compliance defined by equation (2.27) with the gross tooth compliance for the two 
gear teeth 

( yB + yF) / L given by equations (2.15) or (2.18) and (2.25); see equation (2.14). 


2.3.3 Load - Deflection Relationship of a Gear Mesh 

When the total compliance for a single pair of teeth is known for all points of 
contact along the pressure line, computation of the gear-mesh compliance can be carried 
out. Two cases are possible: 

1 . Only one pair of teeth is in contact; then the mesh compliance is equal to the single - 
tooth pair compliance. This is normally the case when contact is near the pitch point. 

2. More than one pair of teeth is in contact. Successive pairs of teeth are spaced along 
the line of action by one normal pitch. Consequently, the compliance curves for 
successive tooth -pairs also can be spaced along the line of action by the same 
amount. 

For simplicity, the real gear - action is modeled as shown in Figure 6 which 
represents engagement and disengagement of the various pairs of teeth. The gear action 
can be represented by the movement of a cam underneath successive pairs of teeth. T his 
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cam is shown in Figure 6 during its passage underneath tooth-pairs n, n-1 and n+1; 
however, after leaving these pairs, it will continue on, bringing other tooth pairs into 
contact in succession. Point P on the cam surface corresponds to pitch-point contact, and 
when any tooth-spring is contacting the cam at P, that tooth-pair is in contact at the pitch- 
point in the real gear system. Tooth-springs in the model are spaced in the gear by an 
interval of one normal pitch, just as the actual tooth-pairs are spaced in the gear system. 
The distance AC on the cam is flat, and is equal to the ideal path of contact. Any tooth- 
pair start engagement at point E and is in full engagement at point A and remain in 
contact until point C is reached, then starts disengagement, untilpoint D is reached. 

The load W is the total load transmitted through the mesh, and Sr is the relative 
motion between mating gears, measured along the pressure line. The relative motion is 

resisted by the tooth-pair stiffiiess, kn , k^, , k^i , etc., depending on the number of teeth 
in contact. 

For example, tooth-pairs n and n+ 1 are in contact at the position shown in Figure 6 , 
the compliance for the mesh at this position is given by the resultant compliance of the 




tooth-springs n and n+1. Since the load L is shared between the two springs, the total 
compliance is the parallel-combination of the individual compliance. 


or 


1 _ 1 1 _ 

c ~ c + r 

total '-'n 


c - 

total 


c c 

c + c 


(2.30) 


(231) 


This result is easily generalized to the combined compliance for n pairs of teeth. 

tic, 


c. = 


r - 1 


nc, Tic, 

r=IL4:=l m=r+\ J 


(2.32) 


where each C refers to the component compliance of a given tooth-pair at a cer tain 
location along the pressure line. 


2.4. Model of Imperfection in Gear Tooth 

To simulate gear feilures, a computer code developed at NASA Lewis Research 
Centerf 1 0 ) was modified to simulate various types of gear tooth imperfection. The 
imperfection of the gear tooth was modeled by using a tooth profile modification 
procedure that simulates changes in gear tooth surface profile. The tooth profile 

modification is represented by their respective cams in Figure 6. The cams representing 

the tooth profile modifications have the form 

(S - S e ) x and Cd (-Sd - S) x respectively (2.33) 

where x is an integer. Typical gear tooth profile modifications were found to be 
represented quite well by simply using a cam form with x - 2; see Figure 6. 
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Hence the imperfection of the tooth can be simulated by engagement and 
disengagement tooth profile modifications Sd and e* , which represent the deviation of the 

profile from the ideal tooth geometry, and are assumed to be a square function of the 
distance along the line of action, e.g. 

£je = C t (Sj -S,) 2 and Sj^CASj-S,) 2 ( 2 . 34 ) 

vdiere Sj reflects the distance along the line of action from the pitch line; Sd and S e , the 
distance along line of action from pitch line to the start of disengagement cam and end of 

engagement cam; C d and C e , disengagement and engagement cam relief, which are 
determined by: 

c e = A e /(S oe -S') 2 and C d = A, /(S^ -S d ) 2 (2.35) 

Where Sod and S<*, represent distances along the line of action from pitch line to end of 
disengagement cam (tooth tip) and start of engagement cam; Ad and A« , maximum 
disengagement and engagement tooth relief for the tooth pair, which is the total or 
combined tooth pair relief at the start and end of an ideal mesh. 

Figure 7 shows how this model in Figure 6 relates to the real gear tooth. 

Suppose the tooth shown in Figure 7 is the driven gear. In Figure 7, the point P is the 
pitch point of the tooth pair; point E is the start point of engagement; point A is the end 
point of engagement; point C is the start point of disengagement; and point D is the end 
point of disengagement. 
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Figure 7. Gear tooth with modified tooth profile 
(The lengths of Soe and S e are measured along the line of action) 


Tooth profile modification can be divided into two parts: one is the engagement 
part, which is the tooth face between points E and A; another is the disengagement part, 
which is the tooth face between points C and D. When two gears are in mesh, modifying 
the tooth tip of one gear is the same as modifying the tooth root of the other gear. 
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Hence the disengagement part can be defined as the driving gear tooth face between point 
A andE. 

To simulate the wear of a tooth face, two parameters needed to be specified. One 
is Psod / Psoe , which is the length of profile modification as a percent of the length of 
total disengagement or engagement measured along the line of action. This parameter 
determines the length of tooth wear. Another one is A<j and A« , which is the maximum 
disengagement and engagement tooth relief for the tooth pair, which is the total or 
combined tooth pair relief at the start and end of an ideal mesh. This parameter 
determines the maximum wear at the tip or the root of the gear tooth. 

2.5. DiscussionsofResults 

As outlined above, there are two parameters which control the simulated tooth 
wear. A parametric study was conducted to evaluate die change in gear-mesh stiffness 
under different wear conditions. The analysis given above, along with its corresponding 
computer code was used to analyze the gear-mesh stiffness change of a gear system under 
different wear conditions. 

For simplicity, a single stage gear mesh was chosen (one input gear and one 
output gear). The pertinent parameters for the system are given in Table 1 . An example 
of the input data for the system is given in Table 2. 
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output gear 


Diametrical Pitch 


8.4667 


8.4667 


Pressure Angle (degree) 22.5 


Number of Teeth 


Face Width (in) 



1.1811 


1.1811 


Input Torque (lb.in) 30 



Input RPM 


Table 2 Gearbox Input Data Set 
Gearbox, Input pinion & gear 


1. 2. 

8.4667 22.5 

0. 

>. 14. 

28. 


>. 1.1811 

1.1811 



30. 

15. 6000. 6000. 1. 

1 28. .0280835 

1 3L .077535 

1 52. .1 

60. 50. 50. .000021 .000 021 

66. 60. 4 0. .00001 

~7Z. !00001 " " 

120 . 0 . 

140. 25. 180. .0528 

150. .01 20. 

167. .0075 

481. .00000000 .00000000 .00000000 
521. .00000 .0000 .00000 

699. 0. ~ 

806. 1. 
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First, Let’s analyze the wear of engagement: 

L Psoe was fixed at 50%, that is, the length of wear was fixed; the tip relief, 

, was chosen as le-5 in, 2e-5 in, and 4e-5 in . The results are shown in 
Figure 8. 

2. was fixed at 2e-5 in, that is, the maximum amount of wear at the tip or 
the root of the gear tooth was fixed; the length of wear, , varied at 
50%, 60% ,70% and 80% . The results are shown in Figure 9, 

Second, Let’s analyze the wear of disengagement part: 

1 • Psod was fixed at 50%, that is, the length of wear was fixed; the tip relief , 

A<i , was chosen as le-5 in, 2e-5 in, and 4e-5 in . The results are shown in 
Figure 10. 

2. Ad was fixed as 2e-5 in, that is the maximum amount of wear at the tip or 
the root of the gear tooth was fixed; the length of wear, P^ , varied as 
50%, 60%, 70%, and 80%. The results were shown in Figure 11. 

Finally, both the wear of the engagement and disengagement were considered: 

1 . Psoe and P^ were fixed at 50%, that is the length of wear was fixed; the tip 
relief, A« and Ad , were chosen as 2e-5 in, 4e-5 in, 6e-5 in, and 8e-5 in. 

The results are shown in Figure 12. 

2. A« and Ad were fixed at 2e-5 in, that is the ma ximum amount of wear at the 
tip or the root of the gear tooth was fixed; the length of wear, P** and P^, 
varied as 50%, 60%, 70%, and 80%. The results are shown in Figure 13. 
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• tip wear 4e-5 (50%) 
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distance along line of action in terms of time t (s) 


Figure 8. Mesh stiffness due to the tooth surface wear on engagement 
(Different value of tip relief, while die wear length remains fixed) 
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Figure 9. Mesh stiffness due to tooth surface wear on engagement 
(Different value of wear length, while the tip relief remains fixed) 
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Figure 10. Mesh stiffiiess due to tooth surface wear on disengagement 
(Different value of tip relief, while the wear length remains fixed) 
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8G% wear (tip. wear 2e-5 in) 
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diistance along line of action in terms of time t (s) 


Figure 1 1 . Mesh stiffiiess due to tooth surface wear on disengagement 
(Different value of wear length, while the tip relief remains fixed) 
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Figure 12. Mesh stiffness due to tooth surface wear on both engagement and 
disengagement (Different value of tip relief, while the wear length remains fixed) 
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3.5e+6 


3.0e+6 


2.5e+6 


2.0e+6 


1.5e+6 


' has no wear 
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Figure 1 3. Mesh stiffiiess due to tooth suiftce wear on both engagement and 
disengagement (Different value of wear length, while tip relief remains fixed! 
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From Figure 8-13, it is obvious that the gear tooth damage due to surface pitting 
and wear can significantly change the phase of the mesh stiffiiess. The higher degree of 
surface pitting and wear , the more the phase of the mesh stiffiiess will shift. 

2.6. Summary 

A numerical model has been developed to simulate the gear mesh stiffiiess change 
resulting from gear tooth damage due to surface pitting and wear. The work in this 
chapter can be summarized as follows: 

A method has been developed to simulate the tooth surface wear in a gear 
transmission system. The tooth surface wear level can be controlled by adjustment 
of both amplitude and length in the tooth profile modification. 

2. The gear mesh model has been developed to provide the mesh stiffiiess with the 
effect of gear tooth damage due to various degrees of surface pitting and wear. It 
will be shown in the next chapter that these changes in gear mesh stiffiiess can be 
incorporated into a dynamic simulation of the gear transmission system for 
dynamic predictions. 
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CHAPTER 3 

THE EFFECTS OF GEAR TOOTH IMPERFECTION ON THE DYNAMIC 
CHARACTERISTICS OF GEAR TRANSMISSION SYSTEMS 


3.1. Objectives 

In the last two decades, problems arising from excessive gear tooth wear and gear 
tooth surface pitting in gear transmission systems have been of increasing concern for a 
variety of gear users. Although regular visual inspections and preventive maintenance 
can help to reduce the failure rate of gear systems, the cost and down time required make 
such programs inefficient and uneconomical. 

Vibration signature analysis methodologies are being developed to non-intrusively 
examine the health and wear of gear transmission systems. Considerable success has 
been achieved in applying the Wigner-Ville distribution concept (WVD) [14-16] in 
machine health monitoring of gear transmission systems [17,18,24]. However , a 
complete vibration signature database is necessary for the development of an effective 
pattern recognition scheme. In order to populate such databases, the development of an 

accurate analytical procedure to predict vibrations in gear systems due to wear and fatigue 
failure is necessary. 

The objective of this chapter is to develop a comprehensive procedure to simulate 
and analyze the vibration in a gear transmission system with effects of surface pitting and 
wear of the gear teeth under normal operating conditions. To simulate the vibration of 
the transmission system, the equations of motion were established individually for each 
rotor-gear-bearing system. The changes of the mesh stiffoess at one particular tooth 
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location or a number of consecutive teeth due to the effects of surface pitting and wear 
were incorporated into the gear-rotor model for the dynamic simulation [1 1-13]. The 
dynamics of each gear-rotor system were coupled with each other through the gear mesh 
interacting forces. The coupling between the rotors and the casing structure were joined 
through the bearing support forces. The global vibrations of the system were evaluated 
by solving the transient dynamics of each rotor system simultaneously with the vibration 
of the casing. In order to minimize the computational effort, the number of degrees-of - 
fioedom of the system were reduced by using a modal synthesis procedure [1 1,12]. The 
results were evaluated by Wigner-Ville Distribution (WVD) [13-15]. 

3.2. Dynamics of the Gear-Shaft Configuration and the Gearbox System 

The dynamics of the ith individual gear-shaft system can be evaluated through the 
equations of motion for the vibrations of an individual rotor-bearing-gear system as 
shown in Figure 14, given in matrix form, as 

[m\w]+ [K,\w t }= fo(/)}+ fo(/)}+ {F\( 0 } (3.1) 

where [M] and [K s ] are respectively the mass and shaft stiflhess matrices of the rotor, 

{ Wj} is the general displacement vector of the ith rotor in the its local coordinate system, 
and* {F w (t)>, (Fgj (t)}, and {Fui (t)} are the force vectors acting on the ith rotor system 
due to bearing forces, gear mesh interactions, and mass-imbalances respectively. 
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Figure 14. Schematic of a rotor-gear bearing system 


The 


equation of motion of the 
Figure 14. Schematic o 


gearbox with p rotor systems can be expressed as 
f a rotor-gear bearing system. 


Ik p . }+ [k 1W =£ [t„ fc , ( o } 

/=! 


(3-2) 


vdiere [TcJ represents the coordinate transformation between die ith rotor and the 
gearbox. 
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The bearing forces {F bi (t)] for the ith rotor can be evaluated as 

to «) = to KW- to ito !)+ to Uk } - to jto }) (3 .3) 

where [C b i ] and [K^ ] are respectively the damping and stiffness of the bearing, [T ic ] is 
the coordinate transformation matrix for the gearbox with respect to the ith rotor, and Wd 
are the casing displacements at the rotor locations. 

The gear forces generated from the gear mesh interaction[40] can be written as 

fc(o}=to(o}+to(o} 0.4) 

where {Fn (t)} is the vector containing the gear forces and moments resul ting from the 
relative rotation between the two mating gears and {F ti (t)} is the vector containing gear 
forces and moments due to the translational motion between the two gears. 

In order to calculate the transient and steady state dynamics of the system, the 
coupled rotor and casing equations of motion must be solved simultaneously. To 
minim ize the computational effort, the modal transformation [1 1,12] procedure is applied 
to reduce the degrees of freedom of the global equations of motion. U sing m undamped 
mode shapes of the ith rotor system [<Jhi , , <J>a , ...» <|>im ] and me undamped mode 

shapes of the gearbox , <|>c 2 » 4*c3 »••>> <j>cmc ]» the rotor displacement for the ith rotor can 

be written as 

W = Z4fe} (3.5) 

and, similarly, the gearbox displacements as 

W=fc]k} 


(3.6) 
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where { Ai } and (A* } are the modal time functions of the ith rotor and the gearbox 
respectively. Using the expansion in equation (3.5), the equations of motion for the ith 
rotor in equation (3.1) can be written as 

[M y, ]fa }+ fa y ]fa } = fa (/)}+ fa, (/)}+ fa (0} (3.7) 

Premultiplying by [<fc ] T and using the orthogonality conditions of the mode shapes [21], 
the ith rotor equations of motion can be written as 

fa fa } = fa, fa fa* fa fa (3.8) 

where [ A 2 ] is the diagonal matrix of the squares of the natural frequencies of the system. 

For the gearbox system, a similar transformation is carried out and equation(3.2) 
can be written as 

RK[A 2 k} = fc} (3.9) 

For the system of k rotors, equation (3.8) can be repeated k times and solved with 
the casing equation (3.9) simultaneously for the modal accelerations (Ai } and {Ac }. A 
numerical integration scheme is used to integrate the accelerations to obtain velocities 
and displacements at each time step for transient calculations [1 1], 

3.3. Vibration Signature Analysis 
3.3.1 Joint Time-Frequency Analysis 

To examine the vibration signal, a joint time-frequency analysis method was 
chosen. This approach was chosen because of the large amount of information 
represented in the joint time-frequency results which can not be represented separately in 
either the time domain or the frequency domain. The joint time-frequency analysis will 
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provide an instantaneous frequency spectrum of the system at every instant of the 
revolution of the pinion while a Fourier Transform can only provide the average vibration 
spectrum of the signal obtained during one complete revolution. In other words, the 
time-changing spectral density from the joint time-frequency spectra will provide 
information concerning the frequency distribution concentrated at that instant around the 
excited instantaneous frequency which cannot be obtained in a regular vibration 
frequency spectrum. The following is a description of the joint time-frequency analysis 
method. 

To examine the vibration signal in a joint time-frequency domain, the Wigner- 
Ville method [14,15] was used in this research. While the Fast Fourier Transform (FFT) 
technique can provide the spectral contents of the time signal, it cannot distinguish time 
phase change during a complete cycle of operation, hi other words, it assumes that the 
time signals are repeatable for each time data acquisition window without considering the 
effects of any magnitude and phase changes during the sampling period. The Wigner- 
ViUe distribution will provide an interactive relationship between time and frequency 
during the period of the time data window. The comprehensive representation of the 
vibration signal using the WVD method is the primary reason that it was used to compare 
the predicted and experimental vibration results. The WVD (Wigner-Ville Distribution), 
in a Discrete form, can be written as: 

W x ( nT , f) = IT 2 x(nT + iT)x\nT - iT ) • (3.10) 


where 
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W(t,f) - the Wigner-Ville distribution in both the time domain t and frequency domain f. 
x(t) = the time signal 

T = the sampling interval 
L = the length of time data used in the transform 

To allow sampling at the Nyquist rate and eliminate the concentration of energy 
around the frequency origin due to the cross product between negative and positive 

frequency [14,15], the analytic signal was used in evaluating the WVD. The analytic 
signal s(t) is defined as 

s{t) = x(t) + jH[x(t)\ 

where H[x(t)] is the Hilber transform of x(t). However, an alternative approach can be 
used to calculate the analytic signal using the frequency domain definition. The analytic 
signal s(t) can be evaluated by calculating the FFT of the time signal x(t), then setting the 
negative frequency spectrum to zero. The analytic signal can be obtained by evaluating 
the inverse FFT of the spectrum. 

To simplify the computational effort, the WVD can be evaluated using a standard 
FFT algorithm. Adopting the convention that the sampling period is normalized to unity, 
it is necessary only to evaluate the WVD at tim e zero. Hence 

L 

K(°J) = 2£*(z> (312) 

i=-L 

where k(i)= s(i) s* (i). 

In order to avoid a repetition in the time domain WVD, a weighting function [28] 
was added to the time data before the evaluation process. Such a process may decrease 
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the resolution of the distribution, but it will eliminate the repetition of peaks in the time 
domain and the interpretation of the result is substantially easier. 

3.3.2 Frequency Domain Analysis 

The frequency spectrum is found by applying a Discrete Fourier Tr ans form (DFT) 
on the time averaged signal x(t), such that the spectral components are 

z ( i) = r|x(O e xp[^] p. 13 ) 

where 

x(t) - time domain signal 
X(k)= frequency domain signal 
T= sampling time interval 
N= number of data points 

The frequency components were examined in the frequency domain. 

3.4. Discussions of Results 

In order to examine the sensitivity of the system vibrations on gear tooth 
imperfections, the vibrations in gear systems due to various levels of gear wear were 
analyzed. The basic parameters used in this analysis are the magnitude and geometry of 
tooth profile deviation and the number of teeth involved. 

The model of rotor-gear system used in this analysis is shown in Figure 15. The 
number of teeth in the gear model is 28. And the vibration analysis of this system under 
various levels of gear wear are given as follows. 
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First, surface pitting on a single tooth is used with the magnitude of the tip relief 

given in Table 3. The damaged tooth is the 12th tooth from the reference point on the 
gear. 


Table m Magnitude of the Maximum Tip Relief 



easel 

case n 

casein 

1st tooth 

l.e-5 

2.e-5 

4.e-5 

4.e-5 

4.e-5 

4.e-5 

4.e-5 

4.e-5 

4.e-5 

2nd tooth 



0 


2.e-5 

4.e-5 

4.e-5 

4.e-5 

4.e-5 

3rd tooth 









4.e-5 


Figure 16 - 18 show the WVD, the time signal (to the left of the WVD), and the 
frequency spectra (below the WVD) of the vibrations of the corresponding system. In 
Figure 16, a small cross pattern had developed in the WVD around the 154 degrees 
location which is the exact mesh time of the 12th tooth. A small phase shift is also 
detected in the time signal at that location. As the damage of the tooth becomes more 
severe (increase in the magnitude of the profile modification), the cross pattern is more 

obvious and the phase shift is getting more pronounced, as shown in Figure 17 and Figure 
18. 

Secondly, the dynamics involved with damage on two consecutive teeth (12th and 
13th) was simulated. During this simulation the wear on file first tooth (12th) is kept 

constant while the amounts of wear on the second tooth (13th) increased similarly to that 
given in Table 3. 
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Figure 19-21 show the WVD, the time signal, and the frequency spectra of the 
vibrations of a corresponding system with damage on two consecutive teeth. When the 
damage on the tooth is increased, the cross pattern shown in the WVD is not as sharp as 
those with single tooth damaged. However the sideband components at the tooth pass 
frequency increase with increasing damage, as seen in the frequency spectra of Figures 
19-21. In addition, the once per revolution low frequency component increases along 
with its sidebands. As seen in the time domain, the phase shift becomes more 
pronounced als.o 

Thirdly, the dynamics involved in three consecutive teeth (12th, 13th and 14th) 
were examined. In this case, the wear/damage in the first and the second teeth are kept 

constant, the amounts of wear in the third tooth increase similarly to those given in Table 
3. 

Figure 22 - 24 show the WVD, the time signal, and frequency spectra of the 
vibrations of the corresponding system with three consecutive teeth damaged. The 
frequency components near the tooth pass frequency show very small changes from those 
with two consecutive teeth damaged. Figure 19-21. However, the frequency components 
near the shaft frequency has acquired a substantial increase with large sidebands. 

Based on the above discussion, one can generalize die effects of single and 
multiple consecutive tooth damage on the vibrations of a gear transmission system. 

Using the WVD technique, 3-D image pattern recognition, procedures can be developed 
to identify various combinations and levels of tooth damage. 
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Gear 1 



Figure 1 5. Gear transmission model 
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Figure 17. Simulated pinion vibration signature 
(only one tooth is damaged, tip relief is 2.0*10' 5 in). 
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Figure 19. Simulated pinion vibration signature 
(two consecutive teeth are damaged, tip relief is 4.0*l<T 5 in, 1.0*10' 5 in respectively) 
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Figure 20. Simulated pinion vibration signature 
(two consecutive teeth are damaged, tip relief is 4.0*10* 5 in, 2.0*10' 5 in respectively). 
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Figure 21. Simulated pinion vibration signature 
(two consecutive teeth are damaged, tip relief is 4.0*10* 5 in, 4.0*10' 5 in respectively) 
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Figure 22. Simulated pinion vibration signature 
(three consecutive teeth are damaged, tip relief is 4.0*10' 5 in, 4.0*10‘ 5 in, 1.0*10 

respectively). 
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Figure 23. Simulated pinion vibration signature 
(three consecutive teeth are damaged, tip relief is 4.0*10' 5 in, 4.0* 10' 5 in, 2. 

respectively). 
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3.5. Summary 

A numerical procedure has been developed to simulate the vibration in a gear 
transmission system with effects of gear tooth damage due to wear and pitting. The work 
presented in this chapter can be summarized as follows' 

1) A modal synthesis methodology has been used to simulate the dynamics of gear 
transmission systems. The computational efforts has been greatly reduced by modal 
transformation. 

2) The gear mesh model developed to simulate the gear tooth damage due to wear and 

pitting can easily be incorporated into the global transmission system for dynamics 
predictions. 

3) The Wigner-Villle distribution (WVD) method provides a compreh ensi ve 

representation of the vibration signal. It was successfully used to verify the analytical 
model. 

4) The WVD method provides detailed information. Hence, using the time averaging 
technique, frequency spectrum analysis, and the WVD, a signature analysis scheme 
can be developed to examine and characterize the vibration signal of the gear system. 

5) A parametric study of the effects on the vibration signal due to various degrees of 
pitting and wear damage, could provide a comprehensive database for gear fault 
detection and damage estimation research. 


CHAPTER 4 


VIBRATION SIGNATURE ANALYSIS OF 
A FAULTED GEAR TRANSMISSION SYSTEM 


4.1 Objective 

The objective of this chapter is to examine and compare the three different 
approaches to detect gear wear and Mure. The frequency domain analysis is based on the 
spectral display from a Fast Fourier Transform algorithm. The time domain analysis 
includes the time synchronous averaged signal, and the statistical based techniques FMO, 
FM4, NA4*. and NB4* applied to the time averaged signal. The joint time-frequency 
analysis uses the WVD on the vibration data with special windowing techniques. All of the 
analysis methods are applied to the vibration data of a Mure from the spiral bevel gear 
fatigue test rig in the NASA Lewis Research Center. Results from the various methods are 
compared and general conclusions are drawn from the results. 

4.2 Technical Approach 

As discussed in the previous section, three major methodologies: A) the frequency domain 
approach, B) the time domain approach, and C) the joint time-frequency approach are used 
in this study. The following is a description of the three methodologies: 

(A) Freq uency Domain Techniques 

The frequency spectrum is found by applying a discrete FFT on the time averaged signal 
x(t), such that the spectral components are 


58 


X(k) - t£x(0 

-« N (4 

where x(t)= time domain signal, X(k)= frequency domain signal, T- sampling time 
interval, N= number of data points. The frequency components are ex amine d in the 
frequency domain and compared with those obtained at various stages of the fault 
development in the spiral bevel pinion. 


(B) Time Domain Techniques 

Four different time domain techniques for early detection of gear tooth damage are used in 
this study for evaluation and comparison. All of the time domain techniques are applied to 
the vibration signal after it has been time synchronously averaged. These techniques are: 
FMO, FM4, NA4*, and NB4*. These parameters are defined as follows: 

FM0[46] is a course fault detection parameter that compares foe maximum peak to peak 
amplitude to foe sum of foe mesh frequencies and its harmonics 


FMO 



(42) 


where PP= maximum peak to peak amplitude in signal, A(fi>= amplitude at mesh 
frequency and harmonics, n= total number of harmonics in frequency range, and FM4[46] 
is an isolated fault detection parameter, and is given by foe normalized kurtosis, of foe 
resulting difference signal as 


y 


FM4 = N £(4-J)V 


r=J 




f = i 


(4.3) 
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where d(t>- A(t) - R(t), _d 2- mean value of d(t), A(t>= original time synchronous signal, 

R(t)= regular meshing components plus their first order side bands, and N = total number 
of data points in the time signal. 

NA4*[45-48] is a general fault detection parameter with trending capabilities. A 
residual signal is constructed by removing regular meshing components from fire time 
averaged signal. For NA4, the first order sidebands stay in the residual signal and the fourth 
statistical moment of the residual signal is then divided by the averaged variance of the 
residual signal, raised to the second power. The average variance is the mean value of the 
variance of all previous records in the run ensemble. This allows NA4 to compare the 
current gear vibration with the baseline of the system under nominal conditions. NA4 is 
given by the quasi-normalized kurtosis equation shown below: 

NA4*(M) = N V(r ( -r) 4 / j— V 

/=. [ M U 

where r = residual signal, r 3= mean value of residual signal. N= total number of data 
points, in one time record, i= data point number in time record, j= time record number, and 
M=current time record number in run ensemble. 

An enhancement to this parameter is given by NA4*, in which the value of foe 
averaged variance is "locked” when foe instantaneous variance exceeds a pie-determined 
value. [47] This provides NA4 with enhanced trending capabilities, in which foe kurtosis of 
foe current signal is compared to foe variance of foe locked baseline signal under nominal 




(4.4) 


conditions. 
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NB4 is another parameter similar to NA4 that also uses the quasi-normalized 
kwtosis given in equation (4.4). The major difference is that instead of using a residual 
signal, NB4 uses the envelope of the signal banpassed about die mesh frequency. Again, as 
with NA4*, NB4* is an enhancement to the NB4 parameter, in which die value of the 
average variance is "locked" when the instantaneous variance exceeds a pie-determined 
value. The equation for NB4* is given below: 


nw*(M)-n£( S , -?) 4 /{.lf; 

M [ M J, x 


N 

l 

i~\ 




and 


(4.5) 


s(t) = magnitude of {b(t) + i {H[b(t)]}} (4 ^ 

where b(t)= time averaged signal bandpassed filtered about the meshing frequency, 
H[b(t)]=the Hilbert Transform of b(t), N= total number of data points in one time record, i— 
data point number in time record, j= time record number, and M= current time record 
number in run ensemble. 

(C) Joint Time-Frequencv Technique 

To examine the vibration signal in a joint time-frequency domain, the Wignen-Ville 
method[14-18,41,49] is used in this study. While the FFT technique (eq. 4.1) can provide 
the spectral contents of the time signal, it cannot distinguish time phase change durin g a 
complete cycle of operation. In other words, it assumes that the time signals are repeatable 
for each time data acquisition window without considering the effects of any magnitude and 
phase changes during die sampling period. The Wigner-Ville distribution will provide an 



interactive relationship between time and frequency during die period of the time data 
window. The WVD (Wigner-Ville Distribution) can be written as: 
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W(t, t) = f>(t+ 1) x* (t - 

_ao “ 


(4,7) 


To allow sampling at the Nyquist rate and eliminate die concentration of energy around die 
frequency origin due to the cross product between negative and positive frequency, [ 1 4, 1 5] 
the analytic signal was used in evaluating the WVD. The analytic signal s(t) is defined as 


Cl 


s(t) = x(t) + jH[x(t)] 

where 

H[x(t)] = the Hilbert transform of x(t) defined by: 

H[x(t)] = 

X.t K 


(4.8) 


(4.9) 


However, an alternative approach can be used to calculate die analytic signal using the 
frequency domain definition. The analytic signal s(t) can be evaluated by calculating die 
FFT of die time signal x(t), then setting the negative frequency spectrum to zero. The 
analytic signal can be obtained by evaluating the inverse FFT of the spectrum. To simplify 
the computational effort, the WVD can be evaluated using a standard FFT algorithm. 
Adopting the convention that the sampling period is normalized to unity, equation (4.7) can 
be rewritten as 


L 

W x (n, f) = 2£x(n+i) x*(n-i) . e' j4,rfi 

i=-L 


(4.10) 


As for the continuous time case, it is necessary only to evaluate the WVD at time zero. 
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Hence 


Wx(0, f) = 2£k(i) e' j4,rfi 

i=-L 


(4.11) 


where k(i) = s(i) s* (-i). Equation (4. 12) can be evaluated using the discrete FFT algorithm. 

In order to avoid a repetition in the time domain WVD, a weighting functional 8] is added 
to the time data before the evaluation process. Such a process may decrease die resolution 
of the distribution, but it will eliminate the repetition of peaks in die time domain and the 
interpretation of the result will be substantially easier. 


4.3 Description of Experimental Procedure 

The fatigue damage on the test pinion shown in Figures 25 to 32 was obtained using 
the spiral bevel gear fatigue test rig illustrated in figure 33, at the NASA Lewis Research 
Center. The primary purpose of this rig is to study the effects of gear tooth design, gear 
materials, and lubrication types on the fatigue strength of aircraft quality gears.[50] Because 
spiral bevel gears are used extensively in helicopter transmissions to transfer power between 
nonparallel intersecting shafts, the use of this fatigue rig for diagnostic studies is extremely 
practical. Vibration data from an accelerometer mounted on the pinion shaft bearing 
housing was captured using a personal computer with an analog to digital conversion board 
and anti-aliasing filter. The 12-tooth test pinion, and the 36-tooth gear have:, 0.5 141 in 
pitch, 35 degree spiral angle, 1 in. face width, 90 degree shaft angle, and 22.5 degree 
pressure angle. The pinion transmits 720 hp at nominal speed of 14,400 rpm. The test rig 
was started and stopped several times for gear damage inspection. The test was ended at 
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17.79 operational hours when a broken portion of a tooth was found visually during one of 
the shutdowns. 

4.4. Discussions of Results 

A series of pictures showing the deterioration of the pinion teeth at various stages of 
the test are illustrated in figures 25 to 32. In figure 1 the initiation of a small pit on one of 
the pinion teeth during the first shutdown, at about five and a half hour into the test is 
shown. As the test progressed, the rig was shut down seven more times to examine the 
severity of the pitting and its relationship with the corresponding vibrations. Figures 26 to 
28 show the increase of the damaged area at fire pinion tooth as the elapsed time increased 
to 6.55, 8.55, and 10.03 hr respectively. Note that in figure 29, as the elapsed time increased 
to 12 hr, the damage of the pinion tooth increased to 75 percent of the tooth surface. At this 
stage, pitting also initiated on the adjacent tooth and continued to grow as the time increased 
to 14,53 hr, figure 30. At 16.16 hr, the damage has grown to three adjacent teeth as shown 
in figure 3 1 . The test was terminated when a breakage is detected on one of fire three 
heavily pitted teeth at 17.79 hr, as shown in figure 32. 

Figure 34 depicts the running speed of the test rig during various stages of die 
experiment. Note that there is some fluctuations present in the running speed after each 
shutdown, with a magnitude of approximately 6 percent about the nominal pinion speed of 
14,400 rpm. There is a sharp change in speed at approximately 8.75 hr. These variations in 
speed create a substantial effect on the vibration signal, which is amplified in die NA4 (fig. 
37), NB4 (fig. 38), and the WVD (fig. 45) analysis. 
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Figure 25. Photograph of Figure 26. Photograph of Figure 27. Photograph of 

pinion damage at 5.50 hr. pinion damage at 6.55 hr. pinion damage at 8.55 hr. 




Figure 28. Photograph of 

Figure 29. Photograph of 

Figure 30. Photograph of 



pinion damage at 10.03 hr. 

pinion damage at 12.03 hr. 

pinion damage at 14.53 hr. 


— 
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Figure 31. Photograph of pinion damage at 
16.16 hr. 


Figure 32. Photograph of pinion damage at 
17.79 hr. 



Accelerometer 
mounting location 
(measuring vertical 
accelerations of 
bearing housing) 



>< 


Figure 33. Spiral bevel gear rig at NASA Lewis Research Center. 




Figure 35 shows the results of the FMO analysis. As seen from the figure, FMO shows 
only moderate changes as the damage starts and progresses. It does not provide any 
indication as the damages extends to the adjacent teeth, resulting in pinion tooth fracture. 
The majority of the variations in the FMO parameter are most probably due to the speed 
changes experienced during the test 

Results from the FM4 parameter, as seen in figure 36, shows a possible reaction as die 
pitting started to occur, however, it does not provide any coherent indication of the severity 
of the pitting as the damage increased. In addition, it does not provide any information to 
distinguish the pitting of a single or multiple teeth. 

Results from NA4 and NA4* are illustrated in figure 37. It is obvious from the figure 
that both NA4 and NA4* provide a very good indication of the pitting development on the 
pinion tooth. The magnitude of the parameter increases to a nondimensional value of 7 after 
shutdown #2 at 6.55 hr, and further to a value of 17 when die pitting covers 75 percent of 
the tooth surface at 8.5 hr. As expected, the "locked" denominator inNA4* provides a more 
robust indication as the pitting progresses.[45»47] Again, both parameters are very sensitive 
to the speed variations, especially after shutdown #3 and #6. 

The NB4 and NB4* parameters, as shown in figure 38, show a very similar trend to 
those of NA4 and NA4*, with a more robust indication to the severity of the damage. 
However, both NA4 and NB4 did not provide any type of indication as the damage spread 
to other teeth, and finally as tooth fracture occurred. 
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Figure 34. Rig operating speed at various run times. 



Figure 35. FMO parameter at various run times. 



Figure 36. FM4 parameter at various run times. 
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Figure 39 shows a WVD (and corresponding intensity scale) and below the WVD the 
frequency spectrum (from FFT analysis) of a uniform sine wave signal shown at the left 
side of the figure. Note that only the frequency component of the input frequency is detected 
by both the frequency domain analysis and the joint time-frequency distribution. The WVD 
does not exhibit any changes during the 1 cycle (0 to 360 degrees) rotation of the shaft 
When a short term amplitude and phase change is added to the system, as shown in figure 
40, the frequency spectrum remains virtually unchanged. The WVD shows a dramatic 
change of die energy distribution pattern at the location where the change occurs. The 
lighter shades of the distribution display indicates a smaller vibration amplitude, which is 
shown by the time signal at the left side of the figure. Such an effect could be possibly 
caused by a chipped or cracked tooth. Figure 41 shows the effects of a short term amplitude 
increase in the time signal to simulate vibrations caused by gear tooth surface damage. Note 
that the frequency spectrum remains the same showing only the component of the exciting 
frequency while the WVD again provides a good indication of die amplitude increase by the 
widening of the shaded area to a diamond shape at the correspon ding "damaged" tooth 
location. Figure 42 shows the effects of a time decaying short term amplitude and phase 
change signal. The WVD shows a half diamond shape of shaded area, similar to that of 
figure 40, at the location where the amplitude and phase changes are presented. As seen in 
figure 42, the frequency spectrum gives very little indication of the signal change. 
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Figure 41. Example of WVD on short term amplitude increase in a signal. 
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Figure 42. Example of WVD on short term amplitude and phase change in a signal. 
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Figures 43 to 50 show the WVD and frequency spectra of the spiral bevel pinion 
vibration at various stages of damage, corresponding to die photographs given in figures 25 
to 32. Figure 43 shows the occurrence of the initial pitting at around 200 degrees from the 
triggering point of the gear, at 5.5 hr, of running time. At 6.55 hr, the pitting on die tooth 
surface progressed to a more noticeable stage, as shown in figure 26, the WVD pattern, 
figure 44, begins to adopt those of a short term amplitude and phase change as illustrated in 
figure 42. At this stage, a Cross pattern appeared in the joint time frequency domain 
(WVD) as the damaged gear tooth produced a change in the frequency components other 
than the mesh frequency. Due to the speed increase at 8.55 hr, the overall WVD amplitude 
increases substantially as shown in figure 45. At the running time of 10.03 hr, the corre- 
sponding WVD in figure 46 shows the initiation of a cross pattern as surface pitting in the 
damaged tooth becomes more pronounced which can also be evidenced from the 
photograph of the pinion gear in figure 28. This phenomenon is also evident in the 
frequency spectrum with the existence of sideband components. At 12.03 hr, when die 
pitting on the pinion tooth has extended to about 75 percent of the tooth surface as shown in 
figure 29, the WVD pattern exhibits a solid cross pattern extending over the mesh frequency 
and several of its adjacent sidebands. The high concentration in the WVD energy and die 
initiation of a second cross pattern at 14.53 hr, figure 48, shows the advancement of the 
pitting process on second tooth, as illustrated in figure 30. This is further confirmed by die 
large amplitude of sideband component(above mesh) in the frequency spectrum. At 16.16 
hr, as seen in figure 49, the WVD pattern changes, showing more advanced damage pattern 
similar to the multiples of the decay of a single short term amplitude increase and phase 
change demonstrated in figure 42. Such phenomenon is due to the pitting of three 
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consecutive teeth in the pinion gear as shown in the picture given in figure 3 1 . The 
frequency spectrum in figure 49 shows a substantial amplitudes increase in the sideband 
components. The discontinuity of the WVD at the mesh frequency, shown in figure 50, 
similar to die example shown in figure 16 due to the short term phase change, is probably 
die result of the instantaneous phase change caused by the fractured tooth, as illustrated in 
figure 32. The two cross patterns in the WVD is very distinct as the effects of pitting at the 
teeth adjacent to the fracture tooth become more pronounced. Note, also that, as given in 
figure 50, the amplitude of the sideband fiequencies(above and below die mesh frequency) 
increase substantially. 

4.5 Conclusions 

Based on the results of the application of the various aforementioned methods, the 
Mowing conclusions can be made: 

1) The FMO parameter shows only moderate changes as the damage starts and progresses. 
It also fails to indicate the fracture of the pinion tooth. 

2) The FM4 parameter shows a possible reaction to the start of the pitting process, 
however, no coherent indication is provided for die growth and severity of the pit. 

3) The NA4* and NB4* parameters show good reactions to die initial pitting damage and 
very nice indications for the growth and severity of the pitting damage. However their 
indications for the tooth fracture is somewhat unclear. 

4) The WVD provides vital information concerning both the severity and the location of 
the pitting process in the gear system. 
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5) The fracture of the gear tooth and its exact location can be pinpointed nsing the WVD 
technique. However a machine vibration signature database is required to interpret the 
resulting WVD. 

(5) The occurrence and the severity of gear tooth failure can be reliably detected using a 
combination of the time averaging, the frequency analysis, and die WVD te chni que 
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CHAPTER 5 

ANALYSIS OF THE EFFECTS OF SURFACE PITTING AND WEAR ON THE 
VIBRATIONS OF A GEAR TRANSMISSION SYSTEM 

5.1 Objective 

Vibration signature analysis methodologies are being developed to non-intrusively 
examine the health and wear of gear transmission systems. Using spectral analysis, the amplitude 
of the frequency spectrum of the measured vibration signal is calculated and displayed in a 
continuous manner. However, the spectral analysis technique is difficult to apply in a highly 
complex system where the large number of spectral lines often makes it difficult to detect 
significant changes in the spectrum. Another methodology is the joint time-frequency approach 
which applies the Wigner-Ville distribution (WVD) [14-16] on the time vibration signal of the 
system. Unlike the regular Fourier transform process, the WVD provides an instantaneous 
frequency spectrum of the system at any instant throughout the sampling period (while FFT 
provides a averaged frequency spectrum of the total sampling period). The spectral density of the 
fundamental exciting frequency and its sidebands change as the shaft rotates through a complete 
revolution. Some success has been achieved in applying the WVD concept in the health 
monitoring of gear transmission systems[16-18,26]. However, a complete vibration signature 
database is needed for development of an effective pattern recognition scheme. In order to 
populate such databases, the development of an accurate analytical procedure to predict 
vibrations in gear systems due to wear and fatigue failure is necessary. 

The objective of this paper is to develop a comprehensive procedure to simulate and 
analyze the vibration in a gear transmission system with effects of surface pitting and wear of the 
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gear teeth under normal operating conditions. The effects of changes in magnitude and ph»s* of 
the mesh stiffness at one particular tooth or a number of consecutive teeth were evaluated in 
order to simulate the effects of surface pitting and wear. The effects of these localized changes in 
the gear mesh were incorporated into each gear-rotor model for the dynamic 
simulation!! 1,51 ,52] . The dynamics of each gear-rotor system were coupled with each other 
through the gear mesh interacting forces. The coupling between the rotors and the casing 
structure were generated through the bearing support forces. The global vibrations of the system 
were evaluated by solving the transient dynamics of each rotor system simultaneously with the 
vibration of the casing. In order to minimize the computational effort, the number of 
degrees -of-freedom of the system were reduced by using a modal synthesis procedure[ll,51]. 

The global transient dynamics of the overall transmission system were calculated in the modal 
coordinates where the modal accelerations were integrated numerically to obtain the velocities 
and displacements at each time step. An FFT procedure was used to transform the transient 
vibrations into the frequency domain for signature analysis. In addition, the Wigner-Ville 
distribution[14-18,26,53] was also used to examine the gear vibrations in the joint 
time-frequency domain for vibration pattern recognition. Experimental vibration results obtained 
from a gear fatigue test rig at NASA Lewis Research Center[48] were used to verify the 
analytical procedure. 

5.2. Solution Procedures 

5.2.1 Dynamics of the Gear-Shaft Configuration and the Gearbox System 
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The dynamics of the ith individual gear-shaft system can be evaluated through the 
equations of motion for the vibrations of a individual rotor-bearing-gear system as shown in 
Figure 51(11,51], given in matrix form, as 

[M] { Wj + [KJ{ Wj = { Fbi(t)} + { Fgi(t)} + { FuA )} (5.1) 

where [M] and [KJ are respectively the mass and shaft stiffness matrices of the rotor, {WJ is the 
general displacement vector of the ith rotor in the its local coordinate system, and, {F^t)}, 
{Fgj(t)}, and (F^t)} are respectively the force vectors acting on the ith rotor system due to 
bearing forces, gear mesh interactions, and mass-imbalances. In this model, the dynamics 
between the gearbox and the rotor are coupled through the bearing forces, which are evaluated 
by the relative motion between the rotor and the gearbox. The interactions between the each 
individual rotor are coupled through the gear forces generated by the relative motion of the two 
mating gears at the mesh point. 

The equations of motion of the gearbox with p rotor systems can be expressed 
as 

[ Me] fwj + [ Kc] {We} = £ [ Ta] {Fbt (t)} (5.2) 

where [T ci ] represents the coordinate transformation between the ith rotor and the gearbox. 

5.2.2 Evaluation of Bearing Forces 

The bearing forces {F bi (t)} for the ith rotor can be evaluated as 
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{ F», (t)} = [ Cj ({ Wj - [ Tic] { Wj) + [ K bi ] ({Wj - [ Tic] { Wet }) (5.3) 

where [CJ and [KJ are respectively the damping and stiffness of the bearing, [TJ is the 
coordinate transformation matrix for the gearbox with respect to the ith rotor, and W d are the 
casing displacements at the rotor locations. 

5.2.3 Evaluation of Gear Forces 

The gear forces generated from the gear mesh interactional can be written as 

<F„(t)} = { Fri(t )} + {Fti(t)} (5.4) 

where (F n (t)} is the vector containing the gear forces and moments resulting from the relative 
rotation between the two mating gears and (F ti (t)} is the vector containing gear forces and 
moments due to the translational motion between tire two gears, the forces and the torsional 
moments due to relative rotation between the kth location of the ith rotor and the 1th loation of 
the jth rotor, can be respectively expressed[54,10] as 


(Fnk(t)} = [ Dgi]{ Kgij} (( Rik 6ik) - (~Rji Oji)) 


(5.5) 


{ Mrik (0} = Ri [ Dgi] { Kgij ((R ik Ovc) - (Rj, 0J,)) (5.6) 

where ^ and 9* are the radius and angular displacements in die ith-rotor gear localized gear 
rotational coordinates, Rj, 7 and $ Jt 8 are radius and angular displacements of the 1th location at 

the jth rotor in the ith-rotor localized gear rotational coordinates, [K,J is the gear mesh stiffness 
matrix between the ith and jth rotors, and [DJ is a diagonalized matrix transforming foe 



localized gear coordinates into the ith rotor coordinate system. In addition, the translational 
forces can be represented by 
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{ F uk (t)} = f D#] {Kgijf [ Dgi /' ([ T,j] { Wj,} - { W ik }) (5.7) 

where [Ty] is the transformation matrix between the ith and jth rotor coordinates. For a rotor 
with one gear mesh, the gear force vector {F g ,(t)} defined by equation (5.4) will have non-zero 
elements only at the gear location, i.e., the kth station. As the gear rotates, the stiffness of the 
gear mesh changes as the gear mesh progresses from single to multiple tooth contact. This 
cycling effect of the gear mesh stiffness is the main source for the mesh frequency excitation in 
the system. 

5.2.4 Modal Synthesis Procedure 

In order to calculate the transient and steady state dynamics of the system, the coupled 
rotor and casing equations of motion must be solved simultaneously. To minimiz e the 
computational effort, the modal transformation[ll,51] procedure will be applied to reduce the 
degrees of freedom of the global equations of motion. Using m undamped mode shapes of the ith 
rotor system [fa, fa, fa. .... fa] and m,. undamped mode shapes of the gearbox [fa, fa, fa, ..., 
fad, the rotor displacement for the ith rotor can be written as 
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(WJ = 2 A, {<t>J 

j-i 

(5.8) 

or 


(WJ = [<f>J{Aj 

and , similarly, the gearbox displacements as 


{Wc} [ <p c ] { Ac} JQj 

where {Aj} and {Aj} are the modal time functions of the ith rotor and the gearbox respectively. 
Using the expansion in equation (8), the equations of motion for the ith rotor in equation (5. 1) 
can be written as 


W [ <f>J { Aj + [ K ,] [ <f>J (AJ = { F « (t)} + { F# (t)} + {Fui (t)} (5.11) 

— Premultiplying by [<fc] T and using the orthogonality conditions of the mode shapes[l 1 ], the ith 

— - rotor equations of motion can be written as 

§P§ 

^ + = (FJ + ff gi) + CfJ ( 5 . 12 ) 

— where [a 2 ] is the diagonal matrix of the squares of die natural frequencies of the system, and 


(Fj = { Fbi(t)} 


(5.13) 
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fFj = [<f>J { Fgi(t)} 


(5.14) 


fFuJ ~ [ 4>i f {Fui(t)j (5.15) 

For the gearbox system, a similar transformation is carried out as equation (2) can be written as 


fAj + [J] {Aj = (Fj (5.16) 

For a system of k rotors, equation (5.12) can be repeated k times and solved with the casing 
equation (5.16) simultaneously for the modal accelerations {AJ and {Ac}. A numerical 
integration scheme is used to integrate the accelerations to obtain velocities and displacements at 
each time step for transient calculations [11]. 

5.3. Signature Analysis of Vibration Signal 
5.3.1 Frequency Domain Analysis 

The frequency spectrum analysis is used by applying a discrete Fourier Transform on the 
average time signal x(t) such that the spectral components are 

x (k) = t Yj x (0 ^P( ~ ) 

i~0 N 


(5.17) 



nj yin j| ( ^ | J jLLuUjulU 
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where x(t) is the tune averaged of the vibration signal W(t) and T is the sampling interval. The 
frequency components are examined in the frequency domain and compared with those obtained 
at various stages of the fault development in the expe rimenta l gear test rig. 




iSt 



S 3 2 Joint Time-Frequency Analysis : The Wigner-Ville Distribution 

To examine the vibration signal in a joint time-frequency domain, the Wigner-Ville 
method[14-16] is used in this study. The Wigner-Ville distribution will provide an inter-domain 
relationship between time and frequency during the period of the time data window. The WVD 
(Wigner-Ville Distribution) can be written as: 


WV(t,j) = { x(t + L) / (t . L) e -r-«fT dr 

-on ^ * 


or in a discrete form as 

WVx (nT.J) = 2T J] x(nT + iT) x '(nT - iT). ? (5.1 9) 

/=*-i 

where WV(t, f) is the Wigner-Ville distribution in both the time domain t and the frequency 
domain f. To allow sampling at the Nyquist rate and eliminate the concentration of energy 
around the frequency origin due to the cross product between negative and positive 
frequeney[14-16], the analytic signal was used in evaluating the WVD. The analytic signal s(t) 
is defined as 


s(t) = x(t) + jH[x(t) ] 

Where H[x(t)] is the Hilbert transform of x(t) defined by : 


(5.20) 
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Hfx(t)] = - J d4 

* i t-z 


(5.21) 


However, an alternative approach can be used to calculate the analytic signal by the frequency 
domain definition. The analytic signal s(t) can be evaluated by calculating the FFT of the time 
signal x(t), then setting the negative frequency spectrum to zero. The analyticsignal can then be 
obtained by evaluating the inverse FFT of the spectrum. To simplify the computational effort, 
the WVD can be evaluated using a standard FFT algorithm. Adopting the convention that the 
sampling period is normalized to unity, Equation (5.19) can be rewritten as 

WV X (n,j) = 2 ^ x(n + i) x (n-i). e ' i4nfi (5.22) 

i=-L 

As for the continuous time case, it is necessary only to evaluate the WVD at time zero. Hence 

WVx (0.J) = 2 2 k(i) e* 4 * (5.23) 

i=-Z, 

where k(i) = s(i) s* (-i). Equation (23) can be evaluated using the discrete FFT algorithm. 

In order to avoid a repetition in the time domain WVD, a weighting function[26] is added 
to the time data before the evaluation process. Such process may decrease the resolution of the 
distribution, but it will eliminate the repetition of peaks in the time domain, and, thus the 
interpretation of the result will be substantially easier. 


5.4. Description of Experimental Study 

The experiment was performed on the spiral bevel gear fatigue test rig [48], as illustrated 
in Figure 52, at the NASA Lewis Research Center. The primary purpose of this rig is to study 
the effects of gear tooth design, gear materials, and lubricants on the fatigue strength of aircraft 
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quality gears. Because spiral bevel gears are used extensively in helicopter transmissions to 
transfer power between nonparallel intersecting shafts, the use of this fatigue rig for diagnostic 
studies is practical. Vibration data from an accelerometer mounted on the pinion shaft bearing 
housing was captured by an analog to digital conversion board. The 12-tooth test pinion, and the 
36-tooth gear have a 35 degree spiral angle, a 1 in. face with, a 90 degree shaft angle, and 22.5 
degree pressure angle. The pinion transmits 720 hpata nominal speed of 14,400 rpm. The test 
rig was stopped several times during the test for gear damage inspection. The test was concluded 
at 17.8 operational hours when a breaken tooth was detected visually during one of the 
shutdowns. 

Pictures of tooth damage on the pinion at various stages in the test are shown in Figure 
53. At the first rig shut-down, at about 5.5 hours into the test, a small pit was observed on one 
of the teeth on the test pinion, as illustrated in Figure 53A. The test was stopped again at 
approximately 12 hours and the pitted area spread to cover approximately 75% of the face of the 
pinion tooth, as seen in Figure 53B. In addition, pitting started to appear on the adjacent teeth. 
Figure 53B shows the pinion at the end of the test, at 17.8 hours. It was found that one of the 
three heavily pitted pinion teeth had experienced a tooth breakage, losing one third of the tooth, 
as shown in the figure. 

5.5. Discussions of Results 

To study the effects of gear tooth pitting and wear on the dynamics of the rotor system, 
the numerical simulation procedure described above was used to model the vibrations of the 
pinion gear in the test rig. During the experimental study, vertical direction vibration signals 
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Figure 51. Schematic of the rotor-gear bearing 
system. 



Figure 52. Spiral bevel gear rig at NASA Lewis 
Research Center. 





Figure 54. Stiffness changes in the gear mesh model. 
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from the pinion gear are time synchronously averaged for spectral analysis and analysis using the 
joint time-frequency distribution(WVD). In order to perform an accurate comparison, the 
averaged time signal from the vertical vibration of the pinion gear is also generated using the 
numerical model. During these simulations, approximate gear mesh stiffness models are 

developed to simulated the effects of wear and pitting of the pinion tooth on the dynamics of the 
system. 

As it has been established, the changes due to gear tooth wear or failure can be 
represented by the amplitude and phase changes in vibration, which, in turn, can be represented 
by magnitude and phase changes in mesh stiffhess[17,18]. To demonstrate the effects of mesh 
stiffness change on gear vibration, the variation of the mesh stiffness model used for this study is 
given in Fig 54. The "undamaged" configuration of the mesh stiffness is given by 0 degree phase 
change (Fig. 54A), and 0% amplitude reduction(Fig. 54B). During the wear and pitting process, 
two types of stiffness changes are examined, i.e., the phase changes, shown in Figure 54A, and 
the amplitude changes, shown in Figure 54B. Figures 55 and 56 show the time, frequency, and 
joint time-frequency analysis(WVD) of the pinion gear vibration signals with the approximated 
changes in gear mesh stiffness. 

Figure 55 shows the effects of mesh stiffness phase changes in the WVD representation 
of the predicted vibration signal. As seen in Figure 55, a phase change in the mesh stiffness at 
the 6th tooth of the 12-teeth pinion resulted in a temporary increase of amplitude and phase of 
the pinion vibration time signal during the 6th tooth pass location. As the phase shift in the mesh 
stiffness increases, from 1.5 degrees to 4.5 degrees, the changes in amplitude and phase in the 
vibration signal become more pronounced. In the frequency spectra, this change in mesh stiffness 
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will result in the increase of the amplitude in the sideband frequencies. However, as discussed 
earlier, although the frequency spectrum provides good indications of the existence of the 
non-synehronous components, it can not distinguish the time locations of their occurrences. The 
joint time-frequency analysis using WVD, as indicated by the scale of shades given in Figure 55, 
shows the existence of various frequency components as the pinion rotates through a complete 
revolution of 360 degrees. Note that, in this case, the WVD shows a continuous excitation of the 
mesh frequency (12 x rotational speed) throughout the complete 360-degree revolution while 

subsynchronous components of 8 times, 4 times, and 1 times rotational speed are occurring at the 
6th and 7th tooth pass locations. 

Figure 56 shows the effects of reduction in mesh stiffness at the 6th gear location. With 
the reduction of mesh stiffness, a substantial change in the vibration at the 6th tooth pass location 
(150 - 180 degrees) is observed. Note that in the case of 50% stiffness reduction, the time 
vibration amplitude at the 6th tooth pass location(at approximately the 180 degree location for the 
12 teeth pinion) almost vanishes and a much larger amplitude at the 7th tooth pass location is 
generated. In addition, the vibration amplitudes of the 8th and 9th tooth pass locations are 
reduced. These reductions in vibration amplitudes at mesh frequency resulted in a much higher 
sub-synchronous components in the frequency spectrum as shown in Figure 56C. The WVD 
shows a distinct type of cross pattern at the intersection of the mesh frequency and the 6th tooth 

pass location with a continuous mesh frequency component throughout the complete pinion 
revolution. 

Figure 57 shows the pinion gear vibration signature analysis of the experimental tinv> 
signal acquired at A) 12 hours when one tooth is severely pitted (Figure 53B), and B) 17.8 




Figure 59 - Numerically simulated pinion vibration due to damage on pinion teeth due to wear and pitting 
(a) Single tooth, (b) Three teeth. 
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hours, when three consecutive teeth are pitted and one has a tooth fracture (Figure 53C). To 
numerically simulate these phenomena, two gear mesh stiffness models, as shown in Figure 58A 
and 58B, which include a combination of phase shift and amplitude change, are introduced. 
Figure 58A represents the gear mesh stiffness for a heavily pitted tooth during a single tooth 
pass. The mesh stiffness is simulated by a 50% loss in stiffness at approximately first 20% of 
contacting period. Figure 58B represent the mesh stiffness, for three tooth pass period, consisting 
of one broken tooth with two heavily pitted teeth at the adjacent sides. Note that the stiffness of 
the middle (broken) tooth is simulated by a 50% loss of stiffness at the first 40% of its contacting 
period while the other adjacent(pitted) teeth are simulated by the stiffness reduction similar to 
that of the single tooth case as shown in Fig. 58A. Additional frictional effects are also added 
into the model to simulate the roughness of the tooth surface due to pitting. The simulated 
vibration signature of the pinion gear is given in Figure 59. Comparing Figures 59A and 57a, for 
the single tooth damage case at 12 hours, one may notice the similarities between both the 
frequency spectra and the WVD display. Some of the unevenness in the experimental time signal 
is mainly due to the modulation of frequencies due to other excitations in the test rig which are 
not numerically modeled. For the tooth break-off case at 17.8 hours. Figures 57B and 59B, both 
the numerical and the experimental WVD display a large cross pattern at the 6th tooth pass 
location due to tooth break-off. However, some discrepancies have been detected between the 
experimental and the numerical time signal at the 4th and 5th tooth pass locations. The 
experimental time signal consists of some higher frequency, smaller amplitude vibration 
modulation, which are not being numerically simulated. This additional modulated signal 
resulted in the excitation of the 14 times rotational speed component, as shown in the frequency 



93 


spectrum in Figure 57B, and, also, in turn, is responsible for the small differences created in the 
WVD. 

5.6. Summary and Conclusions 

A numerical procedure has been developed to simulate the dynamics of gear transmission 
systems with the effects of gear tooth damage due to wear and pitting. The work presented in this 
paper can be summarized as follows: 

1) A modal synthesis methodology has been developed to simulate the dynamics of gear 
transmission systems. While the computational efforts has been greatly reduced by modal 
transformation, the numerical results generated maintains good accuracy. 

2) Gear tooth damage due to wear and pitting can be simulated by amplitude and phase changes 
m the gear mesh stiffness model. The gear mesh model developed can easily be incorporated into 
the global transmission system for dynamic predictions. 

3) Using the time averaging technique, frequency spectrum analysis, and the Wigner-Ville 
distribution, a signature analysis scheme can be developed to e xamine and characterize the 
vibration signal of the gear system. 

4) A parametric study of the effects on the vibration signal due to various changes in the gear 
mesh stiffness model, simulating various degrees of pitting and wear damage, could provide a 
comprehensive database for gear fault detection and damage estimation research. 
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CHAPTER 6 

QUANTIFICATION OF GEAR TOOTH DAMAGE 
BY OPTIMAL TRACKING OF VIBRATION SIGNATURE 


6.1. Objective 

In the last two decades, with demands for higher operating speeds and greater load 
capacity, premature failures in high-performance turbomachinery have often resulted in 
enormous financial losses and, at times, catastrophic consequences. In aeronautical applications, 
where both weight and efficiency are pushed to their design limits, the prevention and 
management of premature equipment failures is a vital part of the maintenance program. Current 
onboard condition-monitoring systems for gas turbine engines often fail to provide sufficient 
time between warning and failure for safety procedures to be implemented. On the other hand, 
inaccurate interpretation of operating conditions may result in false alarms and unnecessary 
repairs and downtime. The early detection of incipient failure in a mechanical system is of great 
practical importance as it permits scheduled inspections without costly shutdowns and indicates 
the urgency and locations for repair before a system incurs catastrophic failure. 

Some success has been achieved in identifying damage in a gear transmission system by 
using the Wigner-Ville distribution (WVD) technique as described in the previous chapters 
[4,7,9,15]. The approach is to use statistical pattern recognition to match the WVD signature 
patterns of damaged gears with standard patterns stored in a data base. Although the WVD 
technique is useful for determining the type and location of the damage, it is not much help in 
quantifying the level of damage. Damage quantification would logically be the next step in 
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failure prediction; however, no explicit attempts at damage quantification have previously 
appeared in the literature. 

This chapter presents a new technique for processing vibration data to quantify the level 
of damage in a gear transmission system. The technique consists of a nonlinear numerical 
optimization in the form of an “optimal tracking” problem [56,57]. The optimization uses a 
dynamic model of the gear mesh and forms an estimate of the time-vaiying mesh stiffiiess that 
best corresponds to the given set of vibration data. The utility of the technique relies on the 
relationship between the wear or damage of a gear tooth and the change in stiffiiess of the mesh 
during a given tooth pass cycle. An analysis of this relationship demonstrates that the 

perturbation of the stiffiiess function from the nominal profile can be used to quantify the level of 
damage. 

The optimal tracking technique for estimating the perturbation of the mesh stiffiiess was 
tested in two settings. First, it was tested on a set of fictitious data generated by computer 
simulation of a one-degree-of-freedom mechanical system with time- varying stiffiiess. The 
solution of the optimal tracking problem matched very closely the actual stiffiiess profile used in 
the model generating the data. Then, the technique was tested on a set of experimental data from 
a gear test rig, but still assuming the one-degree-of-freedom model. Despite the simplicity of the 
model the stiffiiess profile obtained was shown to be useful in correlating to the level of damage 
of the gear transmission system. 

This chapter is organized as follows: Section 2 presents the system model and formulates 
the optimal tracking problem. Section 3 outlines the numerical solution procedure for the 
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nonlinear optimization. Section 4 presents and interprets the results of the optimization and 
discusses the next steps to be taken in developing a comprehensive failure-prediction procedure. 

6 . 2 . OPTIMAL TRACKING PROBLEM 
6.2.1 System Model 

The system considered in this study consisted of a small pinion in mesh with a larger 
gear. A simple model of this system has the two gear masses connected by a spring and a 
damper. The larger gear is much heavier than the pinion; hence, it is assumed to be rigid, so that 
all relative motion between the two is attributed to the motion of the pinion. Then, the equation 
of motion of the pinion takes the form 

mx + cx+ k(t)x = 0 ( 6 . 1 ) 

where m is the mass of the pinion and k(t) and c are the stiffness and damping of the mesh. The 
mesh stiffness is not constant but is nominally a periodic function of the gear angle, with each 
period corresponding to one tooth pass. The high points on the periodic stiffness function 
correspond to gear angles where two pairs of gear teeth are in contact, and the low points 
correspond to angles where only one pair is in contact. 

It has been found in experiments on gearbox vibrations[7-9] that the gear mesh stiffness 
changes with the wear, pitting, or fracture of the gear teeth. Such changes in the gear mesh 
stiffness inevitably lead to changes in the vibration signatures of the mechanical system. The 
objective of the optimal tracking procedure developed in this study is to reconstruct the true 
stiffness profile for a damaged gear tooth from the measured vibration. That is, the objective is to 
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determine the function k(t ) that would result in the measured vibration according to the system 
model (6.1). 

The true stiffness profile can be expressed as the sum of a constant (time averaged) 
component, a nominal periodic component, and a perturbation resulting from gear wear or 
damage. Accordingly, the system model (6.1) is written as 

+ + [ *«. * k periodic - k ptturi>(0 ] * = 0 , (6.2) 

or 

X + —X + n 2 x = u(t)x , (6.3) 

m 

where fP = k aV Q/m and u(t) represents the total time- varying component of the stiffness divided 
by the pinion mass. By defining the variables 

x l =x,x 2 =x, (6.4) 

the system model can be written in the state- variable form 

*i = -& 2 x 2 + u(t)x 2 

m (6.5) 

x 2 = x, 

with the given initial conditions 

*o> ^2 (^0 ) = X 0 . (6.6) 

6.2.2 Optimization Problem 

Suppose that a data set corresponding to the vibration of the pinion is collected over the 
interval [ta,tj\. Let the function describing the data set be denoted as x 2 (t), since it corresponds to 
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the modeled variable x 2 (t). The objective is to determine a reasonable time-varying stiffness 
component u(t) for which the model output x 2 ( t ) approximates the measured data x 2 (t ) . 

A diagram depicting the functional objective is shown in Fig. 60. In the figure u{t) is 
depicted as an input to be chosen so that the error e(t) will be small for all time. Note that this 
problem has the form of a tracking problem, where the control input of a system is designed so 
that the system output follows a prescribed reference function. Such a problem may be 
approached by using the standard techniques of optimal control theory[56,57]. In particular, the 
“design” of a suitable function u(t) may be achieved by minimizing the quadratic cost functional 

2 '/ . 

3 ('/)-?!('/)) + 1 +A» j (/)U, (6.7) 

'o 

where /?,, fa, and /3 f are cost-function weighting coefficients. This form of the cost functional 
penalizes the energy in the error between the modeled output and the measured data. It also 

penalizes the use of too large a stiffiiess perturbation function in order to avoid singularity in the 
solution. 

In the optimal tracking problem the system dynamic equations (6.5) are treated as equality 
constraints imposed in the optimization of the cost (6.7). As such, they are appended to the cost 
function by using time-varying Lagrange multipliers X x (/) and X 2 (t). These Lagrange multipliers 
are themselves governed by differential equations called the costate equations. The costate 
equations together with the state equations of the system model form a two-point boundary value 
problem (TPBVP)[56,57j. The TPBVP equations are 

*i = x, -Q 2 x 2 + u(t)x 2 

m 

*2 = *. 


(State equations) 


( 6 . 8 ) 
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(Costate equations) 1 2 m ' 

X 2 = q 2 a } - u (t); i, - - x 2 (t )) 

(Stationarity condition) 0 = X ,x 2 + fi 2 u(t) 

(Endpoint conditions) x, (t 0 ) = x 0 , x 2 (t 0 ) = x 0 

^i(^) = 0, A 2 (/y) = P f(x 2 (t j) — x 2 (tj) ) % 


(6.9) 

( 6 . 10 ) 

( 6 . 11 ) 

( 6 . 12 ) 


The TPBVP (6.8-6.12) represents a set of necessary conditions for u(t) to be the solution of the 


optimal tracking problem. The TPBVP consists of a set of four coupled differential equations 
(6.8-6.9), together with an algebraic relation (6.10), and some endpoint conditions (6.1 1-6.12) at 
both t 0 and tj. Notice that the TPBVP is nonlinear: the unknown function u(t) multiplies other 
independent variables in the differential equations. 


6.3. NUMERICAL SOLUTION PROCEDURE 

The nonlinear TPBVP (6.8-6.12) is solved by an iterative procedure. A complete and 
general derivation of the procedure is given in Sage[57] and Dyer and McReynolds [55], Some of 
the salient points are outlined below for convenience. 


6.3.1 Successive Sweep Method 

Solving the nonlinear TPBVP requires an iterative method. Although several approaches 
are possible, a common and useful one is to begin with an initial guess u°(t) and use it to 
integrate the nonlinear state equations (6.8) forward in time starting from the initial conditions 
(6.11) to determine the nominal state functions jc{*(r) and x 2 (t). Then, starting from the final 
conditions (6.12), integrate the nonlinear costate equations backward in tim e to dete rmin e the 



100 


nominal costate functions and A 2 °(0- The nominal functions *?(/), A?(t), and 

^2 (0 then satisfy all the TPBVP equations except the stationarity condition (6.10). 

The nominal state, costate, and input functions must be iteratively updated, so that they 
will eventually satisfy all the nonlinear TPBVP equations, including the stationarity condition. 
Each update is accomplished by solving a linearized version of the TPBVP. A standard method 
for doing this is known as the sweep method, whereby a linear relationship between the state and 
costate functions is assumed. Then, the linear TPBVP degenerates into a set of ordinary 
differential equations with endpoint conditions at the final time only. These are solved by a 
straightforward numerical integration. In the case of the optimal tracker these ordinary 
differential equations take the form of the coupled Riccati equations 

c x 2 

Pn = 2 ( — Pw ~Pn ) + Pu-^~ 
m p 2 

Pn=^Pn-P 22 ~Pn( ~ & + «(0 ~ V 2 ) + PuPn~ (6.13) 

m Pi Pi 

Pn - -2p,2 ( -n 2 + „(<) ) + p’i -( fi, + -f- ) 

Pi Pi Pi 

with endpoint conditions 

/>..('/) = A 2 ('/) = 0, Pn(tf) - Pf (6.14) 


together with the auxiliary linear equations 


*1 =h \ ( - + P\\^~ ) ~Pu £ ^r (V2 + Pi<0) 

m p 2 p 2 

( -Q 2 +h(0--^V 2 -A 2 : t- ) “ ^(V 2 +A«(0)CP. 2 * 2 +*i) 
Pi Pi Pi 


(6.15) 
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with the endpoint conditions 

h,(t f ) = h 2 (t f ) = 0. (6.16) 

Note that the x and A variables in the differential equations (6.13) and (6.15) represent the given 
nominal functions. (The zero superscripts have been omitted for convenience.) They are simply 
treated as time-varying coefficients in the numerical integration of the differential equations. The 
solutions of equations (6.13) and (6.15) are then used to compute the corrections to the nominal 
state, costate, and input functions. This computation requires yet another numerical integration, 
this time of the linearized state equations 

Ax, = -Ax, ( ~~ + Pu ~7~ ) + Ax 2 ( -a 2 +«(0--^-V2-P,2-J-) - 
m P 2 P 2 P 2 

x 2 X 

~ h 'lt + +>M0) (6.17) 

P 1 Pi 

Ax 2 = Ax, 

with the zero initial conditions 

Ax,(t 0 ) = Ax 2 (t 0 ) = 0. (6.18) 

Finally, the update of the nominal control is computed as 

Au = -^-(^\ X 1 + ^ +X 2 Q?„AX, +/>12 A *2 +A))> (6.19) 

P 2 Pi 

where s is the step size, and the new nominal control is given by 


w' >, (/)= a'(r) + Au'(f). 


( 6 . 20 ) 
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(The superscripts i and z'+l denoting the iteration number have been reinserted in equation 
(6.20).) The procedure is repeated until the nominal functions converge to a solution. 

The real scalar s e [0,1] in equations (6.15), (6.17), and (6.19) is used as a “step size” 
parameter. Using a smaller value of e tends to decrease the magnitudes of the corrections, thereby 
improving the stability of the iterative procedure but slowing the convergence to the solution. 
Using a larger value of s has the opposite effects. 

6.3.2 Numerical Details 

The choices of the cost-function weighting coefficients /?], pi, and /?f are important for 
effective numerical optimization. The parameter P\ defines the penalty on the difference between 
the calculated and reference vibration signals. Since the goal is to minimize the difference 
between the calculated and tracked vibration signals, a large weighting coefficient (5\ should be 
chosen. The parameter Pi defines the penalty on the function u(t). Generically speaking, Pi 
should impose a lighter penalty on u(t) than P\ imposes on the tracking error. Note also that the 
choice of units has an effect on the appropriate relative sizes of P\ and pi. In the examples 
studied the numerical values of u(t) are considerably larger in magnitude than those of a 
reasonable vibration-signal error; therefore, even if equal weighting between error and control 
were desired, pi should be chosen to be considerably smaller than P\ . An inappropriately large 
choice of the parameter P 2 would make the cost function almost unchanged from one iteration to 
the next. Thus, a small constant value was chosen for the parameter Pi. The parameter /^"defines 



103 


the penalty for the error at the final time. If Pf ’\s too small, a large vibration error at the final 
time will result. 

By following these general guidelines the optimization algorithm described in the 
previous section was realized in a computer program. The equations were integrated with a 
seventh-order Runge-Kutta-Fehlberg method. A summary of the programming steps is given 
below : 

0. Set / = 0 and take the initial guess «0(/) for the function u{t) to be zero. 

1. Using the function u\t) from the previous step, integrate the state equations (6.8) forward in 
time. Calculate the resulting cost function X 

2. Integrate the costate equations (6.9) backward in time. 

3. Use the computed state and costate variables as time- varying coefficients in the integration of 
the Riccati equations (6.13) along with the auxiliary equations (6.15) backward in time. 

4. Integrate the linearized state equations (6.17) forward in time. Using the linearized 
stationarity condition (6.19), calculate the correction A u*(t) to the nominal function v}{t) and 
hence the updated function X"l(r). Also, calculate the new cost function XX 

5. Make decisions about the continuation of the optimization procedure and the choice of the 
parameters: 

a. If the difference between the calculated and tracked vibration signals is small, the 
optimization procedure is finished. 

b. If the difference X + 1 - j < 0 is large enough, repeat from step 1 . 

c. If the difference X + 1 — X < 0 is too small, increase the weighting /3\ and repeat from step 


1 . 
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d. IfJ* +1 > J, repeat from step 1 using a smaller value of the step size s. If this is not 
successful, increase the error weighting f3\ and repeat from step 1. 

Some comments should be made on step 5 of the numerical procedure. It was observed 
that for given values of weighting coefficients and the step-size parameter the optimization 
procedure converges to some value of the cost function. In this case the difference between the 
values of the cost functions > +1 - J 1 becomes negligible after some iterations. This means that 
the cost associated with the control u{t) is becoming dominant. Therefore, it makes sense to start 
a new iteration with an increased weight (5\ (i.e., imposing a higher penalty on the vibration 
error). 

6.4. DISCUSSION OF RESULTS 

To demonstrate the optimal tracking procedure described above, two numerical cases 
were used in this study. The first case was a numerical experiment in which the tracker was 
applied to a set of vibration signals generated numerically, assuming a given gear mesh stiffness 
profile. The mesh stiffness profile evaluated by the optimal tracker was compared with the 
original stiffness used in generating the vibration signal. Figure 61(a) shows the comparison 
between the vibration signal generated by a sinusoidal stiffness and that simulated by the optimal 
tracker. As shown in the figure the two vibration signals were very similar; the small difference 
between the two signals is given in Fig. 61(b). Figure 62(a) depicts the original gear mesh 
stiffness used and the stiffness evaluated by using the optimal tracker; the difference between the 
two stiffnesses is given in Fig. 62(b). The excellent agreement between the two stiffnesses in this 
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numerical experiment has confirmed the applicability of the optimal tracking procedure in 
evaluating system stiffness changes from system vibration signals. However, this close 
resemblance between the generated and simulated signals was partly due to the original time 
signals being smooth, continuous, and harmonic without any substantial change in magnitude 
and phase over the gear revolution. To demonstrate the generality and limitation of the developed 
procedure, a set of experimental data taken from a test rig was used in the next case. 

The second case was based on the experimental data obtained from the spiral bevel gear 
test rig shown in Fig. 63. The primary purpose of this rig is to study the effects of gear tooth 
design, gear materials, and lubrication types on the fatigue strength of aircraft-quality gears [45]. 
Because spiral bevel gears are used extensively in helicopter transmissions to transfer power 
between nonparallel intersecting shafts, using this fatigue rig for diagnostic studies is extremely 
practical. Vibration data from an accelerometer mounted on the pinion shaft bearing housing 
were captured by using a personal computer with an analog-to-digital conversion board and an 
anti-aliasing filter. The 12-tooth test pinion and the 36-tooth gear have the following 
characteristics: 0.5141 in pitch, 35° spiral angle, 1-in. face width, 90° shaft angle, and 22.5° 
pressure angle. The pinion transmits 720 hp at a nominal speed of 14 400 rpm. The test rig was 
started and stopped several times for gear damage inspection. The test was ended at 17.72 
operational hours when a broken portion of a tooth was found visually during one of the 
shutdowns. 

Figure 64(a) depicts the gear tooth after 6.5 hr of operation. Note that there is heavy 
surface pitting on one gear tooth with minor pitting on the next tooth. Figure 64(b) shows the 
time domain averaging, the frequency spectrum, and the joint time-frequency analysis using the 
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evaluated mesh stiffness for 
numerical experiment and error 
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Figure 63 - Experimental spiral bevel 
gear test rig. 
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Figure 64. (a) Photograph of damaged gear at 6.55 hr. (b) Gear vibration signal at 6.55 hr, 
Fast Fourier Transform and Wigner-Ville Distribution. 
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Wigner-Ville distribution (WVD) [4,7-9,15] of the accelerometer signal at 6.5 hr[9]. Note that in 
Fig. 64(b) the time signal indicates a large vibratory signal during the engagement of the sixth 
and seventh teeth (damaged teeth), but the frequency spectrum, because of its averaging 
characteristics, shows very little change from the original signal [9]. The WVD begins to show a 
pattern of shifting of the major frequency component (at a mesh frequency of 2880 Hz) around 
the meshing of the sixth and seventh teeth. The WVD pattern in this case is very similar to those 
resulting from a short-term amplitude and phase change of a vibration signal. Although it has 
been established by the authors in some previous publications [7-9] that such damage in the gear 
can be identified by the WVD pattern recognition process, the level of the damage has not been 
addressed. A recent study by the authors has shown that wear and surface pitting of the gear tooth 
usually will result in a phase shift in the stiffness profile, without any significant change in the 
stiffness magnitude. Figure 65 shows the stiffness change in a gear mesh evaluated [10] from 
gear tooth surface profile variations. Note in Fig. 65(b) that increasing surface profile variation 
increases the phase shift of the gear stiffness without changing the magnitude of the stiffness. 

Incorporating this constant gear mesh stiffness as an additional constraint, the optimal 
tracking procedure was applied to the experimental vibration signal (obtained from the bevel gear 
test rig at 6.5 hr as shown in Fig. 64) to evaluate the corresponding gear mesh stiffness. To better 
evaluate the gear mesh stiffness, the time signal was filtered at a mesh frequency of 2880 Hz. 
Figure 66(a) shows the comparison between the unfiltered experimental signal and the optimal 
tracker simulation, and Fig. 66(b) shows the comparison between the filtered experimental signal 
and the tracker-simulated signal. Note that because of the substantial change of magnitude and 
phase of the time signal during the data acquisition period (one revolution of the gear), the 
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(a) Gear tooth 



Figure 66 - Filtered and unfiltered 
experimental and tracker simulated 
vibration signals for spiral bevel 
gear at 6.5 hr. (a) Unfiltered. 

(b) Filtered. 
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Figure 67 - Tracker evaluated mesh 
stiffness for spiral bevel gear at 


6.5 hr. 



110 


accuracy in the simulated vibration is not as good as that in the numerical experiment (Fig. 
61(a)). Figure 67 depicts the gear mesh stiffness evaluated by using the optimal tracker. Note 
that in the evaluated stiffness considerable phase shifts at several gear teeth resulted in the large 
variation in magnitude and phase of the vibration signal. At the location where pitting occurred 
(teeth 6 and 7) the phase shift of the stiffness was more pronounced. By using the results from 
the evaluated mesh stiffness and the correlation of stiffness change with gear wear shown in Fig. 
65(b), the gear damage can be estimated. 

6.5. CONCLUSIONS 

This chapter presents a unified approach to identifying and quantifying damage in a gear 
transmission system. The conclusions from this study are as follows: 

1. The application of the joint time-frequency technique called the Wigner-Ville distribution 
provides the ability to identify the types and locations of the gear damage. 

2. The optimal tracker developed in this chapter provides a very reasonable estimate of the 
stiffness change due to damage, which can be related to the level of gear damage. 

3. For vibratory signals with large changes in magnitude and phase angle the accuracy of the 
simulated signal from the optimal tracker may decrease. 

4. For a more accurate evaluation of system mesh stiffness an optimal tracker for the complete 
dynamic model of the system is needed. 
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CHAPTER 7 

GENERAL CONCLUSIONS 

During this study, a numerical procedure has been developed to simulate the dynamics of 
gear transmission systems with the effects of gear tooth damage due to wear and pitting. The 
numerical procedure developed consists of (a) a numerical model to simulate the gear mesh 
stiffness change resulting from gear tooth damage due to surface pitting and wear, (b) a mo dal 
synthesis methodology to simulate the dynamics of gear transmission systems, (c) a joint time- 
frequency analysis using the Wigner-ViUle distribution (WVD) method to provide a 
comprehensive representation of the vibration signature for fault detection, and (d) the use of an 
optimal tracker to quantify the damage in the gear tooth. The developed numerical procedure was 
verified by comparison with vibration data from a damaged gear obtained at the NASA Lewis 
Research Center. The specific accomplishements in this project can be summarized as follows: 

1 . A method has been developed to simulate the effects of tooth surface wear in a gear 
transmission system. The tooth surface wear level can be controlled by adjustment of both 
amplitude and length in the tooth profile modification. 

2. Using the developed gear tooth damage model, tooth damage due to wear and pitting can be 
simulated by amplitude and phase changes in the gear mesh stiffness model. The gear mesh 
model developed can easily be incorporated into the global transmission system for dynamic 
predictions. 

3. A numerical procedure has been developed to simulate the vibration in a gear transmission 
system with effects of gear tooth damage due to wear and pitting. The modal synthesis 
methodology was used for simulating the dynamics of gear transmission systems. The gear 
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mesh model developed to simulate the gear tooth damage due to wear and pitting are 
incorporated numerically into the global system for dynamics predictions. 

4. The Wigner-Villle distribution (WVD) method used for the joint time-frequency analysis 
provides a comprehensive representation of the vibration signal. It was successfully used to 
verify the numerical model developed. 

5. Based on the results of the application of the various time and frequency analysis techniques, it 
can be concluded that the WVD provides vital information concerning both the severity and the 
location of the pitting process in the gear system. The fracture of the gear tooth and its exact 
location can be pinpointed using the WVD technique. However a machine vibration signature 
database is required to interpret the resulting WVD. 

6. A parametric study of the effects on the vibration signal due to various degrees of pitting and 
wear damage, has provided a comprehensive database for gear fault detection and damage 
estimation research. 

7. A parametric study of the effects on the vibration signal due to various changes in the gear 
mesh stiffness model, simulating various degrees of pitting and wear damage, could provide 
a comprehensive database for gear fault detection and damage estimation research. 

8. The optimal tracker developed in this chapter provides a very reasonable estimate of the 
stiffness change due to damage, which can be related to the level of gear damage. For 
vibratory signals with large changes in magnitude and phase angle the accuracy of the 
simulated signal from the optimal tracker may decrease. For a more accurate evaluation of 
system mesh stiffness an optimal tracker for the complete dynamic model of the system is 


needed. 
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APPENDIX 


LISTING OF COMPUTER PROGRAM 

FOR 

DYNAMIC SIMULATION 



c************ VLADIMIR V. POLYSHCHUK ****************************** c 
C* *********** poly §fcalpha. mechanical .uakron.edu ************ ******q 
: ^****************** JULY 20, 1995 ***************************C 

^****************** Gear Analysis Program ************************c 
C* NEW DEVELOPMENT CYCLE AS A TRANSITION TO PC WINDOWS C++ PROGRAM C 

^* ************************************* ****Q 
^* ***********************************************(-. 

. J* ***********************************************(-. 

c*************** GUANGHUI XU ************************************c 
:*************** JUNE 7, 1996 *********************************** c 
:*************** modified to include DAMAGED GEAR CASE **********0 
~u*** ******* ************************************* *****************£ 
c* ************************************************************** *£ 

^* ***************************************************************£ 

„ PROGRAM gapd 

PARAMETER N50=50 
, , PARAMETER N20=20 

PARAMETER N10=10 

CHARACTER* 3 0 DATA_FILE (N50) , FSTAT (N50) 

CHARACTER* 30 FILE_OUT (N50) 

CHARACTER* 30 NAME 

— CHARACTER* 1 A, B 

INTEGER NROT, NGEARS, NOBOX, NOMODE, ICOUNT 
INTEGER NMODE (N50) 

INTEGER NUNIT (N50) ,NSTAT(N50) 

INTEGER NBEAR(NIO) , LBEAR (N10 , N10) 

INTEGER NCON BEAR (N10 , N10 , N10 , N10) 
v- INTEGER LONG 

INTEGER NUM, NSWITCH(NIO) 

REAL STIF2 (N10,N10,6,6) , STIFA (N10 ,N10,6,6) 

“ REAL DAMP2 (N10 , N10 ,6,6) 

— REAL SHAPE (N50,N50, N50 , 4 ) 

REAL FREQ (N10,N10, N10) 

~:* ********************************* ********************************** 

v OPEN (9,FILE='menu.dat' ,STATUS='OLD') 

“ DO 100 1=1, N10 

^100 NSWITCH(I) = 0 

— DO 10 1=1, N50 
NUNIT (I ) =9 + I 

^10 CONTINUE 

ICOUNT= 0 
858 READ (9,*) 

w READ (9,*) NROT 

READ (9,*) 

« READ (9,*) NGEARS 

w READ (9,*) 

C* ENTER FILE WITH SHAFT DATA *C 
_ READ (9,15) DATA_FILE (ICOUNT+1) 

FSTAT ( ICOUNT+ 1 ) = • OLD ' 

— READ (9,*) 

C* ENTER FILE WITH GEAR DATA *C 

READ (9,15) DATA FILE ( ICOUNT+2 ) 

FSTAT(ICOUNT+2)= t OLD / ‘ 

T* ENTER FILE WITH BOX DATA *C 
READ (9,*) 

READ (9,*) NOBOX 
_ IF ( NOBOX .EQ. 1 ) THEN 

READ (9,*) 

READ (9,*) 

ICOUNT= ICOUNT - 1 



ELSE 

READ (9,*) 

READ (9,15) DATA FILE (ICOUNT+3 ) 

^ FSTAT (ICOUNT+3 ) = T OLD 7 

ENDIF 
READ (9,*) 

READ (9,*) NOMODE 

— READ (9,*) 

READ (9,15) DATA_FILE ( ICOUNT+4 ) 

READ (9,*) 

READ (9,15) DATA_FILE ( ICOUNT+5 ) 

IF ( NOMODE .EQ. 1 ) THEN 
FSTAT ( ICOUNT+4 ) = ' NEW 7 
_ FSTAT (ICOUNT+5) =' NEW' 

ELSE 

FSTAT ( ICOUNT+4 ) = ' OLD 7 
FSTAT ( I COUNT+ 5 ) = 7 OLD 7 
w ENDIF 

DATA_FILE ( ICOUNT+6 ) = 7 input . out 7 
P FSTAT ( ICOUNT+6 )= 'NEW 7 

NCOUNT=I COUNT+ 6 
DO 30 1=1 , N COUNT 

OPEN(NUNIT (I) ,FILE= DATA FILE (I) ,STATUS= FSTAT ( I ) ) 
“"30 CONTINUE 

Li READ (9,*) 

W READ (9,15) FILE_OUT ( 1) 

READ (9,*) 

— READ (9,*) NUM 

« READ (9,*) ( NSWITCH(I) ,1= 1,NUM ) 

WRITE (NUNIT (NCOUNT) , *) ( NSWITCH (I) , 1= 1,NUM ) 

!p*******************************************************C 
S^******* END OF MENU INPUT ****************************0 
C*******************************************************C 
fe DO 32 1=2 , LEN (FILE_0UT (1) ) 

L, IF ( FILE_OUT(l) (1-1:1) .EQ. 7 7 ) THEN 

LONG= 1-2 
GOTO 33 
r - ENDIF 

S32 CONTINUE 

3 3 CONTINUE 

ij DO 35 1=1 , NGEARS 

~ IF( I .LT. 10 ) THEN 

A= CHAR (48+1) 

0 NAME= FILE_OUT ( 1) 1 

“ OPEN (NUNIT (NCOUNT+I ) ,FILE= NAME ( : LONG) // 7 . 7 //A, 

& STATUS= 7 NEW 7 ) 

ELSE 

~ K= INT (1/10) 

*** J= MOD (1, 10) 

A= CHAR(48+K) 

1 i B= CHAR ( 4 8+ J) 

m NAME= FILE_OUT ( 1 ) 

** OPEN (NUNIT (NCOUNT+I ) ,FILE= NAME ( : LONG) // 7 . 7 //A//B, 

& STATUS= 7 NEW 7 ) 

= ENDIF 

^35 CONTINUE 

^15 FORMAT ( IX, 25A ) 



c*********************************************************************c 

CALL LAT (NUNIT, NCOUNT , NROT , NS TAT , NBEAR, LBEAR, STIF2 , 

& DAMP2,STIFA, FREQ, SHAPE, NMODE,NCON_BEAR, NSWITCH) 

************* ********************************************************q 

IF ( NSWITCH (2) .EQ. 1 ) STOP 

CALL MOTION (NUNIT , NCOUNT , NROT , NSTAT , NBEAR, LBEAR, STIF2 , 

- & DAMP2,STIFA, FREQ, SHAPE, NMODE , NCON_BEAR, NSWITCH) 

C********************************************************************* C 
END 

;********************************************************************* c 

~ SUBROUTINE MOTION (NUNIT, NCOUNT, NROT, NSTAT, NBEAR, LBEAR, STIF2 , 

& DAMP2,STI FA, FREQ, SHAPE, NMODE, NCON_BEAR, NSWITCH) 

_ PARAMETER NONE= 1 

PARAMETER UNBAL= 2 
PARAMETER BEAR= 3 
• PARAMETER GEAR= 4 

v- PARAMETER EXTERN= 5 

PARAMETER ROTOR= 1 
PARAMETER CASING= 2 ’ 

— PARAMETER GROUND^ 3 

PARAMETER N50=50 
_ PARAMETER NROT_SIZE=10 

PARAMETER N10=10 
_ PARAMETER N20=20 

PARAMETER N30=30 
PARAMETER N120=120 
PARAMETER N5 00=5 00 
PARAMETER N5=5 
PARAMETER N1 4 0=140 
~ PARAMETER PI= 3.14 

E :* OLD VARIABLES *C 
^ INTEGER NROT 

w REAL SHAPE (N50 , N50 , N50 , 4 ) 

REAL FREQ (N10,N10,N10) 

3 REAL FTEMP (N10 , N10 , N10) 

INTEGER NMODE (N50) 

m INTEGER NUNIT (N50) , NSTAT (N50) 

» INTEGER NCOUNT 

3 INTEGER NBEAR (N10) , LBEAR(N10,N10) 

652 INTEGER NGEAR(NIO) , LGEAR (N10 , N10) 

INTEGER NCON_BEAR (N10 , N10 , N10 , N10) 

INTEGER NCON_GEAR (N10 , N10 , N10 , N10) 

ii INTEGER NSWITCH (N10) 

REAL STIF2 (N10 , N10 , 6 , 6) , STIFA (N10 , N10 , 6 , 6) 

„ REAL DAMP2 (N10,N10, 6, 6) 

® INTEGER NFACT , NPOINTS , NCALC , NCYCLE , NCYCLE_TRAN 

INTEGER NUNBAL(N20) , ISTAT_UN(N20,N20) 

= INTEGER NUNIT_OUT (N20 , N2 0) 

^ REAL FSAMPLE, TCYCLE(N20) 

REAL UNBX (N20 , N20) , UNBY (N20 , N20) , PHIADD (N20 , N20) 

3* NEW VARIABLES *C 
C* GLOBAL VARIABLES *C 
^ INTEGER NODE_TABLE (N20 , 5 , N50) 

INTEGER CLASS_TABLE(N20,N20) 



INTEGER NODE_CON (N50 , N50) 

INTEGER CLASS_NODE (N50) 

INTEGER CLASS_NAME (N50) 

_ INTEGER NODE_LOC(N50) 

REAL SHAPE_NODE (N50 , N20 , 6) 

REAL STIF_BEAR (N20 , 6,6), DAMP_BEAR (N20 , 6,6) 

— REAL STIF_AV (N20 ,6,6) 

REAL XN(N20 , N20 , 6) , VN (N20 , N20 , 6) , AN(N20,N20 , 6) 

REAL XNP(N20,N20, 6) , VNP (N20 , N20 , 6) , ANP(N20,N20, 6) 

REAL FBEAR(N20,N20, 6) 

— REAL FGEAR (N20 , N20 , 6) 

REAL FUNB (N20 , N2 0,6) 

REAL TORQUE (N20) , FEXT(N20,N20, 6) 

w REAL AMPL_PHASE (N50 , 3 ) 

REAL ANGULA (N20 , 2 ) 

REAL TIME, DELTAT 

INTEGER YES_NO , NL 

INTEGER ITYPE, NTYPE , NT YPE_GLOBAL , NTYPE_STIF 
INTEGER NPOINT (N20) 

— INTEGER ITYPE_DAM(N50,N50) 

REAL AGEAR (N20 , N20 , N10) , AGSTIF (N20 , N120) 

REAL STIF_GEAR (N20 , N120) 

“ REAL TMR(N20 , 3,3) 

REAL ACUR(N20) 

“ REAL GEAR_STIF 

W C AMPL= AMPL_PHASE ( INODE , 1 ) 

0 PHASE= AMPL_PHASE ( INODE , 2 ) 

PHIADD= AMPL_PHASE ( INODE , 3 ) 

C******** INITIALIZATION NC0N_GEAR TABLE *****************************0 
DO 5 IROT=l , NROT 

=; DO 5 IGEAR= 1 , NGEAR ( IROT ) 

w DO 5 JROT=l , NROT 

DO 5 JGE AR=1, NGEAR (JROT) 

51 1 NCON_GEAR ( IROT , IGEAR , JROT , JGEAR) =0 

£ 5 CONTINUE 


n * ******************************************************************** *C 
“;************ INPUT GEAR DATA OF THE ROTOR-BEARING SYSTEM *************0 
^**********************************************************************C 
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INODE= 1 
JCHECK=0 

READ (NUNIT ( 2 ) , * ) 

READ (NUNIT (2),*) 

READ (NUNIT (2),*) NCYCLE , NCYCLE_TRAN , NPOINTS , NFACT 
READ (NUNIT (2),*) 

READ ( NUNIT ( 2 ) , * ) NROT 
READ (NUNIT ( 2 ) , * ) 

READ (NUNIT (2),*) IROT 
READ (NUNIT (2),*) 

READ (NUNIT ( 2 ) , *) ANGULA ( IROT , 1 ) , ANGULA ( IROT , 2 ) 

READ (NUNIT (2),*) 


READ (NUNIT (2),*) NUNBAL(IROT) 

DO 500 I=l,NUNBAL(IROT) 

READ (NUNIT (2 ) , *) ISTAT_UN (IROT , I ) , 

& UNBX ( IROT , I) , UNBY ( IROT , I) , PHIADD ( IROT , I ) 


WRITE (NUNIT (NCOUNT) , *) ' IROT ^IROT, 

& 

WRITE (NUNIT (NCOUNT) ,*) UNBX ( IROT , I ) , 


• NUNBAL ' ,NUNBAL(IROT) , 
' UNBX # , ' UNBY ' 

UNBY (IROT, I) 


“500 CONTINUE 


READ (NUNIT ( 2 ) , *) 

_ READ (NUNIT ( 2 ) , * ) NODE_LOC ( INODE) , TORQUE (IROT) 

CLASS_TABLE ( IROT , EXTE&N) = 1 
NODE_TABLE ( IROT , EXTERN , 1 ) = INODE 
CLASS_NODE( INODE )= IROT 

- INODE= INODE + 1 

READ (NUNIT (2),*) 

READ (NUNIT(2) , *) CLASS_TABLE (IROT, NONE) 

- DO 505 INONE= 1, CLASS_TABLE (IROT, NONE) 

READ (NUNIT (2 ) , *) ILOC 

NODE_LOC ( INODE ) = ILOC 
^ NODEJTABLE ( IROT , NONE , INONE ) = INODE 

CLASS_NODE ( INODE ) =IROT 
INODE= INODE + 1 
505 CONTINUE 

JCHECK=JCHECK+ 1 

IF ( J CHECK .LT. NROT ) GOTO 10 

- READ (NUNIT (2 ) , * ) 

READ ( NUNIT ( 2 ) , * ) YES_NO 
IF ( YES_NO .EQ. 0 ) GOTO 510 

w WRITE (NUNIT (NCOUNT) ,*) 

WRITE (NUNIT (NCOUNT) , *) 7 GEAR DATA INPUT 7 

WRITE (NUNIT (NCOUNT) , *) 

JCHECK=0 

515 READ (NUNIT (2),*) 

READ (NUNIT (2 ) , *) IROT 

- READ (NUNIT (2),*) 

READ (NUNIT (2),*) NGEAR(IROT) 

- DO 520 IGEAR= l,NGEAR(IROT) 

“ READ (NUNIT (2),*) LGEAR ( IROT , IGEAR) 

w WRITE (NUNIT (NCOUNT) , *) LGEAR (IROT, IGEAR) 

READ (NUNIT (2 ) , * ) ( AGEAR (IROT, IGEAR, I) , 1=1,5 ) 

M WRITE (NUNIT (NCOUNT) , *) ( AGEAR (IROT, IGEAR, I) , 1=1,5 ) 

^520 CONTINUE 

JCHECK=J CHECK+ 1 

_ IF ( JCHECK .LT. NROT ) GOTO 515 

WRITE (NUNIT (NCOUNT) , *) 7 GEAR CONNECTION TABLE 7 

WRITE (NUNIT (NCOUNT) , *) 7 IROT IGEAR JROT JGEAR 

H & ITYPE CON 7 

w READ (NUNIT (2),*) 

READ (NUNIT (2),*) NL 

m DO 525 1= 1 , NL 

g READ (NUNIT ( 2 ) , * ) IROT, IGEAR, JROT, JGEAR, ITYPE 

NCON_GEAR ( IROT , IGEAR , JROT , JGEAR) = ITYPE 
^ NCON_GEAR (JROT, JGEAR, IROT, IGEAR) = NCON_GEAR( IROT, IGEAR, JROT, JGEAR) 

5^ WRITE (NUNIT (NCOUNT) ,*) IROT, IGEAR, JROT, JGEAR, ITYPE 

CC & NCON GEAR (IROT, IGEAR, JROT, JGEAR) 

■525 CONTINUE 

: READ (NUNIT (2),*) 

w READ (NUNIT (2 ) , * ) 

“ JCHECK=0 

" READ (NUNIT (2) , *) NTYPE_STIF, NTYPE DAM 

530 READ (NUNIT ( 2 ) , * ) 

READ (NUNIT (2),*) ITYPE 

_ READ (NUNIT (2),*) NL, ' AGSTIF ( ITYPE , 1 ) , AGSTIF ( ITYPE , 2 ) 
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“ AGSTIF ( ITYPE , 3 ) =0 . 0 

NPOINT ( ITYPE ) = NL 

READ (NUNIT(2) , *) (AGSTIF (ITYPE , 3+1) , STIF_GEAR (ITYPE , I ), 1=1 , NL) 

“ WRITE (NUNIT (NCOUNT) ,*) 

WRITE (NUNIT (NCOUNT) , *) ' DATA FOR GEAR TYPE ' , ITYPE 

WRITE (NUNIT (NCOUNT) ,532) 

- & (AGSTIF (ITYPE, 3+1) , STIF_GEAR(ITYPE, I) ,I=1,NL) 

532 FORMAT (IX, 2F20.4 ) 

~ JCHECK= J CHECK + 1 

IF ( J CHECK .LT. (NTYPE_STIF+NTYPE_DAM) ) GOTO 530 
NTYPE_GLOBAL=l 

= NTYPE= NTYPE_STIF + NTYPE_GLOBAL 

READ (NUNIT (2),*) 

i_ DO 535 ITYPE=1 , NTYPE 

READ (NUNIT (2 ) , *) 

READ (NUNIT (2 ) , *) ( (TMR( ITYPE, IROW, JCOL) , JC0L=1, 3) , IR0W=1, 3) 

— WRITE (NUNIT (NCOUNT) , *) 

WRITE (NUNIT (NCOUNT) , *) ' DATA FOR COORD. TRANSFORM TYPE ' , ITYPE 

WRITE (NUNIT (NCOUNT) ,533) 

& ( (TMR( ITYPE, IROW, JCOL) ,JC0L=1,3) ,IR0W=1,3) 

”“ 533 FORMAT ( IX, 3F10.4 ) 

535 CONTINUE 

c*******************-*******************c 

C**** START READING GEAR DAMAGE MODEL *C 
; DO 1999 1=1, NTYPE 

-C DO 1999 J=1 , N30 

C NDAM ( I , J ) = I 

: 1999 CONTINUE 

READ (NUNIT (2),*) 

READ (NUNIT (2),*) NTYPE_DAM 
READ (NUNIT (2),*) 

WRITE (NUNIT (NCOUNT) , *) ' DAMAGE MODEL ' 

DO 2200 1=1 , NTYPE_DAM 

READ (NUNIT (2),*) ITYPE, ITER, NEW_ITYPE 
NDAM (ITYPE, ITER )= NEW_ITYPE 

WRITE (NUNIT (NCOUNT) , *) ITYPE, ITER, NDAM (ITYPE, ITER ) 

C 2200 CONTINUE 
J CHECK = 0 

write (66,*) 'before reading of damage' 

READ (NUNIT (2),*) 

READ (NUNIT (2),*) ITYPE, NTEETH 
DO 1000 1=1, NTEETH 

ITYPE_DAM ( ITYPE , I ) = ITYPE 
CONTINUE 

write (66,*) 'before reading of ndam' 

READ (NUNIT ( 2 ) , * ) NDAM 
WRITE (NUNIT (NCOUNT) ,*) 

WRITE (NUNIT (NCOUNT) ,*) 'DAMAGE DATA FOR GEAR TYPE', ITYPE 
WRITE (NUNIT (NCOUNT) , *) ' ITEETH ITYPE_DAM' , ' NDAM ', NDAM 

IF ( NDAM .LE. 0) GO TO 1020 
DO 1010 1=1, NDAM 

READ (NUNIT ( 2 ) , * ) ITEETH, ITYPE_IDAM 
ITYPE_DAM( ITYPE, ITEETH) = ITYPE_IDAM 

WRITE (NUNIT (NCOUNT) , *) ITEETH, ITYPE_DAM( ITYPE, ITEETH) 
CONTINUE 

J CHECK = JCHECK + 1 

IF (JCHECK .LT. NTYPE_STIF) GO TO 1050 


:c 

' 1050 


1000 

cc 


1010 

1020 


lu 





— * END OF GEAR INPUT DATA FROM ABASE . DAT *C 
510 CONTINUE 

********************************************************C 
u********************************************************c 
DO 700 IROT=l , NROT ‘ 

DO 700 IMODE=l , NMODE ( IROT) 

FTEMP ( IROT , IMODE , 1 ) = FREQ ( IROT , IMODE , 1 ) * 0 . 1047195 
FTEMP ( IROT , IMODE , 2 ) = FREQ (IROT, IMODE, 1) *0 . 1047195 
FTEMP ( IROT , IMODE , 3 ) - IfREQ ( IROT , IMODE ,1)*0.1047195 
FTEMP (IROT , IMODE , 4 ) = FREQ ( IROT , IMODE , 1) *0 . 1047195 
_ FTEMP ( IROT , IMODE , 5 ) = FREQ (IROT, IMODE , 3 ) 

FTEMP ( IROT , IMODE , 6 ) = FREQ ( IROT , IMODE , 4 ) 

700 CONTINUE 

— WRITE (NUNIT (NCOUNT) , * ) 

ii DO 710 IROT=l , NROT 

DO 710 IMODE=l, NMODE (IROT) 

DO 710 ICOORD=l , 6 

FREQ ( IROT , IMODE , ICOORD) = FTEMP ( IROT , IMODE , ICOORD) 

WRITE (NUNIT (NCOUNT) , *) ' FREQ ', FREQ (IROT, IMODE, ICOORD) 

__7 10 CONTINUE 

WRITE (NUNIT (NCOUNT) ,*) 

:************************** A*****************************c 
vj* TABLES FOR CLASSES AND NODES *C 
C* NODE NUMBERING *C 

WRITE (NUNIT (NCOUNT) ,*) ' INODE, IROT, NODE_LOC (INODE) ' 

— WRITE ( NUNIT ( NCOUNT ) , * ) 

DO 650 IROT= 1 , NROT 

WRITE (NUNIT (NCOUNT) ,*) ' BEARINGS NUMBERING ' 

r* NODE NUMBERING FOR BEARINGS 

CLASS_TABLE (IROT, BEAR) = NBEAR(IROT) 

— DO 660 IBEAR=1 , NBEAR ( IROT) 

NODE_TABLE (IROT, BEAR, I BEAR )= INODE 
CLASS_NODE (INODE) = IROT 
NODE_LOC( INODE )= LBEAR ( IROT , I BEAR) 

w WRITE (NUNIT (NCOUNT) ,*) INODE, IROT, NODE_LOC (INODE) 

INODE= INODE + 1 
660 CONTINUE 

I j ’’ 

“ WRITE (NUNIT (NCOUNT) , *) ' GEARS NUMBERING ' 

C* NODE NUMBERING FOR GEARS 

CLASS_TABLE (IROT, GEAR) = NGEAR (IROT) 

S DO 665 IGEAR=1, NGEAR (IROT) 

NODEJTABLE ( IROT , GEAR, I GEAR) = INODE 
CLASS_NODE (INODE) = IROT 
NODE_LOC (INODE) = LGEAR ( IROT , IGEAR) 

— WRITE (NUNIT (NCOUNT) , *) INODE, IROT, NO DE_LOC (INODE) 
INODE= INODE + 1 

665 CONTINUE 

“ WRITE (NUNIT (NCOUNT) , *) • UNBALANCE NUMBERING # 

r* NODE NUMBERING FOR UNBALANCE 

CLASS_TABLE ( IROT , UNBAL) = NUNBAL ( IROT) 

_ DO 670 IUNB=1, NUNBAL (IROT) 

NODE_TABLE ( IROT , UNBAL , IUNB) = INODE 

CLASS_NODE( INODE )= IROT 

NODE_LOC( INODE )= ISTAT_UN (IROT, IUNB) 



AMPL_PHASE (INODE, 1)= SQRT ( UNBX (IROT, IUNB) * 

& UNBX (IROT, IUNB) + 

“ & UNBY ( IROT ; IUNB) *UNBY ( IROT, IUNB) ) 

AMPL_PHASE ( INODE , 2 ) = AT AN 2 (UNBY ( IROT , IUNB) , UNBX ( IROT , IUNB) ) 
AMPL_PHASE( INODE, 3)= £HIADD ( IROT , IUNB) 

~ WRITE (NUNIT (NCOUNT) ,*) INODE, IROT, NODE LOC(INODE) 

INODE= INODE +1 
670 CONTINUE 

“650 CONTINUE 

NNODE= INODE- 1 

_* CHECK 

NNODE= INODE- 1 

DO 222 INODE2=l , NNODE ' 

WRITE (NUNIT (NCOUNT) ,223) ( AMPL_PHASE (INODE2 , I) , 1=1 , 3 ) 

*-223 FORMAT ( IX, 3( E15.7,1X ) ) 

222 CONTINUE 

_* CAN BE MODERNIZE TO RESEMBLE CASING FROM ROTOR 

DO 652 IROT=l , NROT 
CLASS_NAME ( IROT) = ROTOR 
L652 CONTINUE • 

WRITE (NUNIT ( NCOUNT) ,* ) 

WRITE (NUNIT (NCOUNT) , *) ' NNODE = ', NNODE 

MODE SHAPES 

DO 620 INODE=l, NNODE 
IROT= CLASS_NODE( INODE) 

— ILOC= NODE_LOC( INODE) 

DO 620 IMODE=l , NMODE ( IROT) 

SHAPE_NODE ( INODE , IMODE , 1 ) = SHAPE ( IROT , IMODE , ILOC , 1 ) 

_ SHAPE_NODE ( INODE , IMODE , 2 ) = SHAPE (IROT, IMODE, ILOC, 2) 

SHAPE_NODE ( INODE , IMODE , 3 ) = SHAPE ( IROT , IMODE , ILOC , 1 ) 
SHAPE_NODE ( INODE , IMODE , 4 ) = SHAPE ( IROT , IMODE , ILOC , 2 ) 

‘ SHAPE_NODE (INODE, IMODE, 5) = SHAPE ( IROT , IMODE , ILOC , 3 ) 

— SHAPE_NODE (INODE, IMODE, 6) = SHAPE (IROT, IMODE, ILOC, 4) 

WRITE (NUNIT (NCOUNT) , *) ' INODE ', INODE 

IZ WRITE (NUNIT (NCOUNT) ,622) (SHAPE_NODE ( INODE, IMODE , I) , 1=1 , 6) 

“"620 CONTINUE 

WRITE (NUNIT (NCOUNT) ,* ) 

_ WRITE (NUNIT (NCOUNT) , *) ' INODE IROT ILOC ' 

DO 630 INODE=l, NNODE 
IROT= CLASS_NODE (INODE) 
is ILOC= NODE_LOC( INODE) 

5= WRITE (NUNIT (NCOUNT) , *) INODE,' ',IROT,' ' , ILOC 

630 CONTINUE 

™; * **************************************************** * ****** c 
X* NODE CONNECTIONS FOR BEARINGS *C 
WRITE (NUNIT (NCOUNT) , * ) 

~ WRITE (NUNIT (NCOUNT) ,*) ' NODE CONNECTIONS ' 

_ WRITE (NUNIT (NCOUNT) , *) ' INODE JNODE ICON ' 

ICON=l 


DO 621 IROT=l , NROT 
* NODE CONNECTIONS FOR BEARINGS 
DO 621 IBEAR=1 , NBEAR ( IROT) 

DO 629 JROT=l , NROT 

IF ( JROT .GT. IROT ) GOTO 629 



DO 629 J BEAR= 1 , NBEAR ( JROT ) 


IF ( NCON_BEAR(IROT,IBEAR, JROT, JBEAR) .EQ. 1 ) THEN 

— INODE= NODE_TABLE (IROT, I BEAR, I BEAR) 

NODE_CON ( INODE, 1)= 1 

JNODE= NODE_TABLE (JROT , JBEAR , JBEAR) 

NODE_CON ( INODE , 2 ) = JNODE 

— NODE_CON ( INODE , 3 ) = ICON 

WRITE (NUNIT (NCOUNT) , *) INODE, JNODE, ICON 
* SYMMETRY OF THE CONNECTIONS 
NODE_CON (JNODE, 1)= 1 

— NODE_CON (JNODE , 2 ) = INODE 
NODE_CON (JNODE , 3 ) = ICON 

WRITE (NUNIT (NCOUNT) , *) JNODE, INODE, ICON 
w DO 624 1=1,6 

DO 624 J=l, 6 

STIF_BEAR ( ICON, I , J) = STIF2 (IROT, IBEAR, I , J) 

DAMP_BEAR ( ICON , I , J) = DAMP2 ( IROT , IBEAR , I , J) 

— STIF_AV (ICON, I , J) = STIFA(IROT, IBEAR, I, J) 

624 CONTINUE 

ICON= ICON + 1 
i - i ENDIF 

*C* DEVELOP CLEAR CONNECTION TABLE, HERE 2= GROUND 

IF ( NCON_BEAR(IROT,IBEAR, JROT, JBEAR) .EQ. 2 ) THEN 
INODE= NODE_TABLE ( IROT , BEAR , IBEAR) 

_ NODE_CON ( INODE, 1)= 1 

"5* GROUND CONNECTION 
JNODE= NNODE+1 

NODE_TABLE ( NROT+ 1 , BEAR , 1 ) = JNODE 

— CLASS_NODE (JNODE) = NROT+1 

CLASS_NAME ( NROT+ 1 ) = GROUND 

NODE_CON ( INODE , 2 ) = JNODE 

— NODE_CON ( INODE , 3 ) = ICON 

WRITE (NUNIT (NCOUNT) , *) INODE, JNODE, ICON 
DO 625 1=1,6 
DO 625 J=l, 6 

STIF_BEAR(ICON, I, J) = STIF2 (IROT, IBEAR, I, J) 

DAMP_BEAR ( ICON , I, J) = DAMP 2 ( IROT , IBEAR , I , J) 

STIF_AV(ICON, I , J) = STIFA (IROT, IBEAR, I, J) 

_625 CONTINUE 

ICON= ICON + 1 
ENDIF 

5*629 CONTINUE 
621 CONTINUE 

NCON= ICON-1 

^ik**********************************************************^ 

C* CHECK 

g WRITE (NUNIT (NCOUNT) , *) ' BEARINGS NCON = ', NCON 

IS WRITE ( NUNIT ( NCOUNT ) , * ) 

WRITE (NUNIT (NCOUNT) ,*) ' STIF BEAR (IROT, I, J) ' 

DO 71 ICON=l , NCON 

WRITE (NUNIT (NCOUNT) ,72) ( (STIF BEAR (ICON, I ,J) , 1=1 , 6 ) , J=1 , 6) 
-71 CONTINUE 

WRITE (NUNIT (NCOUNT) ,*) 

WRITE (NUNIT (NCOUNT) ,*) 

WRITE (NUNIT (NCOUNT) ,*) ' STIF AV(IROT,I,J) ' 

DO 73 ICON=l , NCON 

~ WRITE (NUNIT (NCOUNT) , 72 ) ( (STIF AV (ICON, I , J) , 1=1, 6) , J=l, 6) 

SJ73 CONTINUE 

WRITE (NUNIT (NCOUNT) , *) 

WRITE (NUNIT (NCOUNT) , *) 



WRITE (NUNIT (NCOUNT ) ,*) ' DAMP_BEAR (IROT, I, J) ' 

DO 74 ICON=l , NCON 

WRITE (NUNIT (NCOUNT) ,72) ( (DAMP BEAR (I CON, I, J) , 1=1,6) , J=l, 6) 
_ 74 CONTINUE 1 

WRITE (NUNIT (NCOUNT) ,*) 

i 72 FORMAT ( IX , 6(E10.4,1X) ) 

- J* 


1* ************************ *********************q 

:* CHECK 

— WRITE (NUNIT (NCOUNT) ,*) ' MODE SHAPES FOR ALL NODES ' 

WRITE (NUNIT (NCOUNT) , *) ' NNODE = ', NNODE 

DO 627 INODE=l, NNODE 
DO 627 IMODE=l , NMODE ( 1 ) 

"" WRITE (NUNIT (NCOUNT) ,622) ( SHAPE_NODE ( INODE , IMODE , I ) ,1=1,6) 

622 FORMAT ( IX, 6(F11.5) ) 

: 627 CONTINUE 

— WRITE (NUNIT (NCOUNT) , *) 

^ 'l**************** **************************** **************Q 

i Z * ****************************************** *c 

*-C* NODE CONNECTIONS FOR GEARS IN THE SYSTEM *C 

ICON= NCON + 1 
DO 750 IROT= 1 , NROT 
~ DO 750 IGEAR= 1 , NGEAR (IROT) 

DO 755 JROT= l,NROT 
IF ( JROT .EQ. IROT ) GOTO 755 
_ DO 755 JGEAR= 1 , NGEAR (JROT) 

IF ( NCON_GEAR ( IROT, IGEAR, JROT, JGEAR) .NE. 0 ) THEN 
INODE= NODE_TABLE ( IROT , GEAR , IGEAR) 

C* FIX IT LATER= MAKE MULTIPLE CONNECTIONS EASY 
NODE_CON ( INODE , 1 ) =1 

— JNODE= NODE_TABLE (JROT , GEAR , JGEAR) 

ITYPE= NCON_GEAR ( IROT , IGEAR, JROT , JGEAR) 

w NODE_CON ( INODE , 2 ) = JNODE 

NODE_CON ( INODE , 3 ) = ICON 
“ NODE_CON ( INODE , 4 ) = ITYPE 

NODE_CON ( INODE , 5 ) = ITYPE 
C* SYMMETRY OF THE CONNECTIONS 
NODE_CON (JNODE , 1 ) =1 ’ 

H NODE_CON (JNODE , 2 ) =INODE 

— NODE_CON (JNODE , 3 ) =ICON 
NODE_CON (JNODE , 4 ) = ITYPE 

rj. NODE_CON (JNODE , 5 ) = ITYPE 

WRITE (NUNIT (NCOUNT) , *) * GEARS CONNECTION 7 

w WRITE (NUNIT (NCOUNT) , *) INODE, JNODE, ICON, ITYPE 

ICON= ICON + 1 
“ ENDIF 

755 CONTINUE 
. 750 CONTINUE 
_ NCON2= I CON- 1 -NCON 

WRITE (NUNIT (NCOUNT) ,*) 

WRITE (NUNIT (NCOUNT) , *) ’ THERE ARE / ,NCON2, / GEAR CONNECTIONS ' 

WRITE (NUNIT (NCOUNT) ,*) 

Q* ************************************** £ 

ml* MODAL EXTERNAL TORQUE CALCULATION *C 
******************** * ****************** C 

DO 760 IROT=l , NROT 

INODE= NODE_TABLE(IROT, EXTERN,!) 



DO 760 IMODE=l , NMODE ( IROT) 

FEXT ( IROT , IMODE , 6 ) = TORQUE ( IROT) *SHAPE_NODE ( INODE , IMODE , 6 ) 
760 CONTINUE 

C* ******************************* *C 

- C* FILE OUTPUT NUMBERING *C 

r; K=1 

- DO 24 IROT=l , NROT 

NNONE= CLASS_TABLE ( IROT , NONE ) 

DO 24 INONE= 1 , NNONE 
NUNIT_OUT ( IROT , INONE ) = K 
~ K= K + 1 

24 CONTINUE 

^***************************c 

DO 205 IROT=l , NROT 
TCYCLE (IROT) =60. 0/ANGULA( IROT, 1) 
r 205 CONTINUE 

^ 2 ** ******************** 

TIME=0 . 0 

NCALC= NCYCLE*NPOINTS 
IF ( NCALC .EQ. 0 ) STOP 
w NCYCLE_OUT=NCYCLE-NCYCLE_TRAN 

NCALC OUT= (NCYCLE-NCYCLE_TRAN) *NPOINTS 
NTOTAl> NCYCLE*NPOINTS*NFACT 
G* GOOD CONVERGENCE IS ABOUT DELTAT=10E-6 *C 
“ FSAMPLE= ANGULA (1,1) *NPOINTS/60 . 0 

4000 CONTINUE 

3 DELTAT= TCYCLE ( 1 ) / ( NPOINTS *NFACT) 

~C IF ( DELTAT .GT. IE-6 ) THEN 

C NFACT= NFACT*2 

2 GOTO 4000 

-C ENDIF 

WRITE ( NUNIT ( NCOUNT ) , * j ' DELTAT= ' , DELTAT , ' NFACT= 7 , NFACT 

~ OPEN ( 4 0 , FI LE= ' head . out ' , STATUS= ' NEW ' ) 

WRITE (40,*) NCYCLE_OUT , NCALC_OUT , 

& ANGULA (1,1) , FSAMPLE 

^ OPEN (41, FILE= ' head2 . out ' , STATUS= 7 NEW ' ) 

WRITE (41,*) NCYCLE_OUT , NCALC_OUT , 

& ANGULA (1,1), NPOINTS 

«i,3****** main loop ****************c 

DO 100 ICALC= 1, NCALC 
*2 DO 105 I FACT = 1, NFACT 

— DO 110 IROT=l , NROT 

CC IF ( ICALC .GT. 2 ) STOP 

rCC PRINT * , 7 BEFORE BEARF ' 

S CALL BEARF ( IROT, NODE_TABLE, CLASS_TABLE, NODE_CON, 

& CLASS_NODE, NMODE, 

& SHAPE_NODE , STIF_BEAR, STIF_AV, DAMP_BEAR , 

S & XN, VN, FBEAR ) 

CC PRINT *, 7 BEFORE UNBF' 

CALL UNBF ( IROT, NODE_TABLE , CLASS_TABLE , NMODE, 

& SHAPE_NODE , ANGULA, AMPL_PHASE , TIME, FUNB ) 

CC PRINT * , ' BEFORE GEARF ' 

m A CUR ( IROT) = 360 . 0*TIME/TCYCLE (IROT) 

“ CALL GEARF ( IROT, NODE_TABLE , CLASS_TABLE, NODE_CON, 

& CLASS_NODE , NMODE, 

& SHAPE_NODE, XN, VN, FGEAR, 

— & AGEAR, ACUR, STIF_GEAR, NPOINT, AGSTIF, TMR, 



~ & GEAR_STIF,ITYPE_DAM ) 

C PRINT *, 'BEFORE EQU' 

_ CALL EQU ( IROT, NMODE, FREQ, FUNB, FBEAR, FGEAR, 

& XN, AN ) 

DO 210 IMODE= 1, NMODE (IROT) 

— DO 210 ICOORD=l , 6 

VNP ( IROT , IMODE, ICOORD) = VN ( IROT, IMODE, I COORD) 

XNP ( IROT , IMODE , ICOORD) = XN ( IROT , IMODE , ICOORD) 

_210 CONTINUE 

CALL EULER ( IROT, NMODE, DELTAT, XNP, VNP, AN ) 

_C PRINT *, 'BEFORE BEARF2 ' 

CALL BEARF ( IROT, NODE_TABLE , C LAS S_T ABLE , NODE_CON, 

& CLASS_NODE, NMODE, 

& SHAPE_NODE, STIF_BEAR, STIF_AV, DAMP_BEAR , 

— & XNP, VNP, FBEAR ) 

■•'C PRINT *, 'BEFORE UNBF2 ' 

CALL UNBF ( IROT, NODE_TABLE, CLASS_TABLE, NMODE, 

— & SHAPE_NODE , ANGULA, AMPL_PHASE , TIME+DELTAT, FUNB ) 

C PRINT * , ' BEFORE GEARF2 ' 

ACUR ( IROT ) = 3 60.0* (TIME+DELTAT) /TCYCLE (IROT) 
w CALL GEARF ( IROT, NODE_TABLE , CLASS_TABLE , NODE_CON, 

& CLASS_NODE, NMODE, 

; : & SHAPE_NODE , XNP, VNP, FGEAR, 

w & AGEAR, ACUR, STIF_GEAR, NPOINT, AGSTIF, TMR, 

& GEAR_STIF, ITYPE_DAM ) 

— CALL EQU ( IROT, NMODE, FREQ, FUNB, FBEAR, FGEAR, 

— & FEXT, XNP, ANP ) 

CALL NEWMARK( IROT, NMODE, DELTAT, XN, VN, AN, ANP ) 

DO 200 IMODE=l, NMODE (IROT) 

DO 200 ICOORD=l , 6 
fig AN (IROT, IMODE, ICOORD) =0.0 

g ANP (IROT, IMODE, ICOORD) =0.0 

200 CONTINUE 

110 CONTINUE 

— TIME=TIME + DELTAT 

105 CONTINUE 

= IF ( ICALC .LE. NCYCLE TRAN*NPOINTS ) GOTO 100 

CALL OUTP ( NUNIT , NCOUNT , NUNIT_OUT , 

S & NROT, NMODE, CLASS_TABLE , 

= & NODE_TABLE , SHAPE_NODE, XN, VN ) 

r. WRITE (6 6,*) GEAR_STIF 

~100 CONTINUE 

c***********************************************************************c 

= RETURN 

— END 

C* ************************************************************ **********c 
m SUBROUTINE BEARF ( IROT, NODE_TABLE , CLASS_TABLE , NODE_CON, 

_ & CLASS_NODE, NMODE, 

& SHAPE_NODE, STIF_BEAR, STIF_AV, DAMP_BEAR, XN, VN, FBEAR ) 


PARAMETER NONE= 1 



PARAMETER UNBAL= 2 
PARAMETER BEAR= 3 
PARAMETER GEAR= 4 

PARAMETER N20=20 
PARAMETER N50=50 

-J* GLOBAL VARIABLES *C 

INTEGER NODE_TABLE (N2 0 , 5 , N50 ) 

INTEGER CLASS_TABLE (N20 , N20) 

INTEGER NODE_CON (N50 , N50) 

" INTEGER CLASS NODE(N50) 

INTEGER NMODEXN50) 

w REAL SHAPE_NODE (N50 , N20 ,6) 

REAL STIF_BEAR(N20 , 6 , 6) , STIF_AV(N20 , 6 , 6) 
REAL DAMP_BEAR (N2 0,6,6) 

REAL XN(N20 , N20 , 6) , VN(N20,N20, 6) 

~ REAL FBEAR (N20 , N20 , 6) , FB(6) 

• C* LOCAL VARIABLES *C 

INTEGER I BEAR, NBEAR, IMODE, I COORD 

REAL XI (6) , VI (6) , X2(6), V2(6) 

;J* INITIALIZE VARIABLES *C 

DO 17 IMODE=l , NMODE ( IROT) 

DO 17 ICOORD=l,6 

w FBEAR ( IROT , IMODE , I COORD ) = 0.0 

17 CONTINUE 

; J* START CALCULATIONS *C 

C WRITE (50,*) ' IN BEAR ' 

gj NBEAR= CLASS_TABLE ( IROT , BEAR) 

DO 100 IBEAR= 1, NBEAR 


SS DO 10 ICOORD=l , 6 

XI (ICOORD) = 0.0 
VI (ICOORD) = 0.0 
10 CONTINUE 

Vir' 

INODE= NODE_TABLE (IROT, BEAR, IBEAR) 

y WRITE (50,*) ' INODE INODE 

DO 110 IMODE=l, NMODE (IROT) 

“ DO 110 ICOORD= 1,6 

^ XI ( ICOORD) = XI (ICOORD) + 

& XN(IROT, IMODE, ICOORD) *SHAPE_NODE (INODE, IMODE, ICOORD) 
„ VI (ICOORD) = VI (ICOORD) + 

_ & VN(IROT, IMODE, ICOORD) *SHAPE_NODE( INODE, IMODE, ICOORD) 

-110 CONTINUE 


— NCON= NODE_CON ( INODE, 1) 

T WRITE (50,*) ' NCON ' , NCON 

; J* LOOP FOR ALL CONNECTIONS *C 


DO 120 ICON= 2, 2*NC0N+1, 2 
JNODE= NODE_CON ( INODE , ICON) 



JCON= NODE_CON (INODE , ICON+1) 
JROT= CLASS_NODE ( JNODE) 


WRITE (50,*) ' IROT INODE JNODE JCON JROT ' 

C WRITE (50,*) IROT, INODE, JNODE, JCON, JROT 

:*? FIND CLASS OF THE NODE *C 

— DO 20 ICOORD=l , 6 
X2 ( ICOORD) = 0.0 
V2 (ICOORD) = 0.0 

___ FB( ICOORD) =0.0 

20 CONTINUE 

DO 130 IMODE= 1 , NMODE ( JROT ) 
w DO 130 ICOORD= 1,6 

X2 (ICOORD) = X2 (ICOORD) + 

& XN (JROT, IMODE, ICOORD) *SHAPE_NODE (JNODE, IMODE, ICOORD) 

V2 (ICOORD) = V2 (ICOORD) + 

— & VN( JROT, IMODE, ICOORD) *SHAPE_NODE (JNODE , IMODE , ICOORD) 

130 CONTINUE 

DO 150 ICOORD= 1,6 
w SUMX=0 . 0 

SUMV=0 . 0 
~ SUMA=0 . 0 

L DO 160 JCOORD= 1,6 

SUMX= SUMX + 

& STIF_BEAR( JCON, ICOORD, JCOORD) *( Xl(JCOORD) - X2 (JCOORD) ) 
SUMV= SUMV + 

— & DAMP_BEAR(JCON, I COORD, JCOORD) *( VI (JCOORD) - V2 (JCOORD) ) 

SUMA= SUMA + • 

& STIF_AV(JCON, ICOORD, JCOORD) *X1( JCOORD) 

160 CONTINUE 

— FB ( ICOORD) = SUMV + SUMX - SUMA 
150 CONTINUE 

WRITE ( 50 , * ) ' FB 1234' 

"C WRITE (50 , *) FB ( 1) , FB(2), FB(3), FB(4) 

* DO 170 IMODE=l, NMODE (IROT) 

FBEAR (IROT, IMODE, 1)= FBEAR (IROT, IMODE, 1) + 

& SHAPE_NODE (INODE, IMODE, 1) *FB(1) + 

& SHAPE_NODE (INODE, IMODE, 2) *FB(2) 

~ FBEAR ( IROT , IMODE , 2 ) = FBEAR ( IROT , IMODE , 1 ) 

li FBEAR ( IROT , IMODE , 3 ) = FBEAR ( IROT , IMODE , 3 ) + 
n & SHAPE_NODE( INODE, IMODE, 3) *FB(3) + 

& SHAPE_NODE (INODE, IMODE, 4) *FB(4) 

FBEAR ( IROT , IMODE , 4 ) = FBEAR ( IROT , IMODE , 3 ) 


W FBEAR ( IROT , IMODE , 5 ) = FBEAR ( IROT , IMODE , 5 ) + 

& SHAPE_NODE (INODE, IMODE, 5) *FB(5) 

FBEAR ( IROT , IMODE , 6 ) = FBEAR (IROT, IMODE, 6) + 

— & SHAPE_NODE (INODE, IMODE, 6) *FB(6) 

IP 170 CONTINUE 

E.. .. 

^120 CONTINUE 
riOO CONTINUE 

c***********************************************************************c 

U RETURN 

g END 


C* ********************************************************************* *c 
^***********************************************************************0 
SUBROUTINE UNBF ( IROT, NODE_TABLE , CLASS_TABLE, 

& NMODE , 

& SHAPE_NODE , ANGULA, AMPL_PHASE, TIME, FUNB ) 

PARAMETER NONE= 1 
PARAMETER UNBAL= 2 
PARAMETER BEAR= 3 
PARAMETER GEAR= 4 

PARAMETER N20=20 
PARAMETER N50=50 


PARAMETER ID_X=1 
PARAMETER ID_Y=3 
PARAMETER PI=3. 1415926 


~C* GLOBAL VARIABLES *C 

INTEGER NODE_TABLE (N20 , 5 , N50) 
INTEGER CLASS TABLE (N20 , N20) 

_ INTEGER NMODEXN50) 

REAL SHAPE_NODE (N50 , N20 , 6) 

REAL FUNB (N20 , N20 , 6 ) 

— REAL TIME 

REAL AMPL_PHASE(N50, 3) 

REAL ANGULA (N2 0,2) 

C* LOCAL VARIABLES *C 

INTEGER IUNB, NUNB, IMODE, ICOORD 

’ w REAL AMPL, PHASE, PHIADD 

REAL ARG 

REAL SPEED, ANGACL 
„ REAL OMEGA 

„ <7* *********************** *c 


SPEED= ANGULA (IROT, 1) * 

ANGACL= ANGULA (IROT, 2) 

w OMEGA=PI*SPEED/30 . 0 

DO 10 IMODE=l, NMODE (IROT) 

L- DO 10 ICOORD* 1,6 

FUNB ( IROT , IMODE , ICOORD) =0 . 0 
FUNB (IROT , IMODE , ICOORD) =0 . 0 
^ 0 CONTINUE 

NUNB= CLAS S_TABLE ( IROT , UNBAL ) 

,_C PRINT *, NUNB 

DO 20 IUNB=1 , NUNB 

E INODE= NODE_TABLE( IROT, UNBAL, IUNB) 

.~C PRINT *, INODE 

t:.:. 

— AMPL= AMPL_PHASE (INODE, 1) 

PHASE* AMPL_PHASE ( INODE , 2 ) 

— PHIADD* AMPL_PHASE ( INODE , 3 ) 

ARG* OMEGA*TIME + PHASE + PHIADD 



cc 


PRINT *, AMPL, PHASE, ARG, OMEGA, TIME 


_ DO 30 IM0DE=1 , NMODE ( IROT) 

FUNB ( IROT , IMODE , ID_X) = 

& AMPL*SHAPE NODE (INODE, IMODE, ID_X) * 

- & ( ANGACL*SIN7ARG) + OMEGA*OMEGA*COS (ARG) ) 

FUNB ( IROT , IMODE , ID_Y) = 

& AMPL*SHAPE NODE (INODE, IMODE, ID_Y) * 

& ( ANGACL*C0S7ARG) + OMEGA*OMEGA*SIN (ARG) ) 

-CC PRINT *, SHAPE_NODE (INODE, IMODE, ID X), FUNB (IROT, IMODE , ID X) 

30 CONTINUE “ " 

20 CONTINUE 

n*** *************** ******************************************************£ 
RETURN 
_ END 

C************************************************'k*'k'ki'it'k'k*ic'k*ic'k-k-k*1t**'k-k1'-k-k*C 

'********-k*******-k*'k*'k'k***'k****'kic****-k***-k**'k-k-k****'k'k*********'k'k*'k*i'ic'kC 
^********************** SUBROUTINE EQUATION *************************c 

C* ******************************************************************* * c 

SUBROUTINE EQU( IROT, NMODE, FREQ, FUNB, FBEAR, FGEAR, 

_ & FEXT, XN, AN ) 

PARAMETER N10= 10 
PARAMETER N20= 20 
L PARAMETER N50=50 

V* GLOBAL VARIABLES *C 

- INTEGER IROT 
INTEGER NMODE (N50) 

REAL FREQ (N10 , N10 , N10) 

~ REAL FUNB(N20,N20, 6) ' 

REAL FBEAR (N20,N20, 6) 

M REAL FGEAR (N20,N20, 6) 

□ REAL FEXT (N20 , N20 , 6 ) 

REAL XN (N20 , N20 , 6) 

REAL AN (N20 , N20 , 6) 

\-4 

w* LOCAL VARIABLES *C 

INTEGER ICOORD, IMODE 

_************* THE EQUATIONS OF MOTION ****************************c 

C WRITE (50,*) ' IN EQU • 

_ DO 100 IMODE=l , NMODE ( IROf) 

AN (IROT, IMODE, 1)= FUNB (IROT, IMODE, 1) - 
& FBEAR (IROT, IMODE, 1) - 

= & FREQ(IROT, IMODE, 1) **2 *XN( IROT, IMODE, 1) + 

& FGEAR (IROT, IMODE, 1) 

^ AN (IROT, IMODE, 2) = AN ( IROT, IMODE , 1) 

AN (IROT, IMODE, 3)= FUNB ( IROT , IMODE , 3 ) - 
& FBEAR ( IROT , IMODE , 3 ) - 

_ & FREQ (IROT, IMODE, 3 )**2*XN (IROT, IMODE, 3) + 

& FGEAR ( IROT , IMODE , 3 ) 

AN ( IROT , IMODE , 4 ) = AN (IROT, IMODE, 3) 



AN ( IROT , IMODE , 5 ) = -FBEAR(IROT, IMODE, 5) - 
& FREQ ( IROT, IMODE, 5) **2*XN (IROT, IMODE, 5) + 

_ & FGEAR ( IROT , IMODE , 5 ) 

AN (IROT, IMODE, 6) = -FBEAR( IROT, IMODE, 6) - 
& FREQ(IROT, IMODE, 6) **2*XN (IROT, IMODE, 6) + 

- & FGEAR ( IROT , IMODE , 6 ) + 

& FEXT ( IROT , IMODE , 6 ) 

WRITE (50,*) ' IROT IMODE ICOORD FUNB FBEAR ' 

-C WRITE (50,*) IROT, IMODE 

C WRITE (50 , *) FUNB(IROT, IMODE, 1) , -FBEAR (IROT, IMODE, 1) 

_100 CONTINUE 

r**************************************************************c 
RETURN 
_ END 

c******************************************************************c 

SUBROUTINE EULER ( IROT, NMODE, DELTAT, XN, VN, AN ) 

— PARAMETER N20=20 
PARAMETER N50=50 

_ INTEGER NMODE (N50) 

REAL H, DELTAT 

_ REAL XN(N20,N20, 6) , VN (N20 , N20 , 6) , AN (N20 , N20, 6) 

REAL XNP(N20,N20, 6) , VNP (N20 , N20 , 6) 

INTEGER IMODE, ICOORD 

C* ************** THE EULER METHOD **************************** *C 
H=DELTAT 

DO 10 IMODE = 1, NMODE (IROT) 

M DO 20 ICOORD=l , 6 

u VNP (IROT, IMODE, ICOORD j= VN (IROT, IMODE, ICOORD) + 

& ’ H*AN (IROT, IMODE, ICOORD) 

XNP(IROT, IMODE, ICOORD) = XN (IROT, IMODE, ICOORD) + 

11 & 0 . 5*H* ( VNP(IROT, IMODE, ICOORD) + VN( IROT, IMODE, ICOORD) ) 

^20 CONTINUE 

10 CONTINUE 

ri DO 30 IMODE= 1, NMODE (IROT) 

DO 30 ICOORD=l , 6 

VN ( IROT , IMODE , ICOORD) = VNP ( IROT , IMODE , ICOORD) 

XN (IROT, IMODE, ICOORD) = XNP ( IROT , IMODE , ICOORD) 

*_3 0 CONTINUE 

RETURN 
L END 

^******************************************************************0 

c* ************************************************************ *c 

SUBROUTINE NEWMARK( IROT, NMODE, DELTAT, XN, VN, AN, ANP ) 

^ PARAMETER N20=20 

PARAMETER N50=50 


INTEGER NMODE (N50) 
REAL H, DELTAT 



REAL XN(N20,N20, 6) , VN (N20 , N20 , 6 ) , AN(N20,N20, 6) 
REAL ANP (N20 , N20 , 6 ) 


— INTEGER IMODE, ICOORD 
REAL BETA 

BETA=0 .167 

C******* NEWMARK-BETA METHOD ************************************** c 
H=DELTAT 

DO 20 IMODE= 1 , NMODE ( IROT ) 

DO 20 ICOORD=l , 6 

— XN ( IROT, IMODE, ICOORD )=XN( IROT, IMODE, ICOORD) + 

& VN(IROT, IMODE, ICOORD) *H + 

& (0.5 - BETA ) *AN (IROT, IMODE, ICOORD) *H*H + 

& BETA*ANP ( IROT , IMODE , ICOORD) *H*H 

— VN(IROT, IMODE, ICOORD) =VN(IROT, IMODE, ICOORD) + 

& 0 . 5*H* ( AN (IROT, IMODE, ICOORD) +ANP (IROT, IMODE, ICOORD) ) 

20 CONTINUE 

RETURN 
II END 

1 _j** ************* ************************* **********************q 

C* ***************************************************** *c 
r.*************** SUBROUTINE OUTPUT *********************c 
,"?*** OUTPUT AND POST-PROCESSING OF THE RESULTS ******c 
^ SUBROUTINE OUTP(NUNIT,NCOUNT,NUNIT OUT, 

& NROT, NMODE, CLASS_TABLE , 

& NODE_TABLE , SHAPE_NODE, XN,VN ) 

— PARAMETER NONE= 1 
PARAMETER UNBAL= 2 
PARAMETER BEAR= 3 

__ PARAMETER GEAR= 4 

PARAMETER N10=10 
i- PARAMETER N20=20 

\Z PARAMETER N50=50 

PARAMETER SC=10000 

INTEGER NUNIT (N50 ) 

— INTEGER NUNIT_OUT (N2 0 , N2 0) 

INTEGER NMODE (N50) 

INTEGER JSTAT, IUNIT 

il INTEGER NROT, IROT, INODE, IMODE 

INTEGER NODE_TABLE (N20 , 5 , N50) 

INTEGER CLASS_TABLE (N20 , N20) 

_ REAL SHAPE_NODE (N50 , N20 , 6) 

REAL XN(N20,N20, 6) , VN (N20 , N20 , 6) 

■ REAL XD(6) , VD (6) 

— IROT=0 
IUNIT=0 

1 DO 10 IROT=l , NROT 

NNONE= CLASS_TABLE (IROT, NONE) 

DO 10 INONE=l , NNONE 

CO 30 ICOORD=l , 6 
XD( ICOORD) =0.0 
VD ( ICOORD) =0 . 0 


30 


CONTINUE 


INODE= NODEJTABLE ( IROT , NONE , INONE ) 

DO 20 IMODE=l , NMODE ( IROT) 

DO 20 ICOORD=l , 6 

XD ( ICOORD) = XD(ICOORD) + 

& XN(IROT, IMODE, ICOORD) *SHAPE NODE (INODE, IMODE, ICOORD) 
VD (ICOORD) = VD( ICOORD) + - 1 

& VN ( IROT , IMODE , ICOORD) *SHAPE NODE ( INODE , IMODE , ICOORD) 

20 CONTINUE ~ ' ’ 


IUNIT= NUNIT_OUT ( IROT , INONE ) 

WRITE (NUNIT (NCOUNT+IUNIT) , * ) XD ( 3 ) 

100 FORMAT ( IX, 2(G15.9,1X) ) 


10 CONTINUE 


RETURN 

END 


;**************************************************************** 

vj******************************************************^*^^^^^^^^ 

C* CODE FOR GEAR FORCE SUBROUTINE 

'.**************************************************************** 

_ SUBROUTINE GEARF(IROT; NODE TABLE, CLASS TABLE, NODE CON, 

& CLASS_NODE , NMODE, “ “ ' 

& SHAPE_NODE , XN, VN, FGEAR, 

& AGEAR, ACUR, STIF GEAR, NPOINT, AGSTIF, TMR, GEAR STIF, 

- & ITYPE_DAM ) ' - ' 


*******0 

*******(2 


PARAMETER NONE= 1 
PARAMETER UNBAL= 2 
PARAMETER BEAR= 3 
PARAMETER GEAR= 4 

PARAMETER N10=10 
PARAMETER N20=20 
PARAMETER N50=50 
PARAMETER N120=120 

PARAMETER PI=3. 1415926 


* GLOBAL VARIABLES *C ] 

— INTEGER NODEJTABLE (N2 6, 5,N50) 

INTEGER CLASS_TABLE (N20 , N20) 

INTEGER NODE_CON (N50 , N50) 

INTEGER CLASS_NODE (N50 ) 

INTEGER NMODE (N50) 

INTEGER NPOINT (N20) 

M INTEGER ITYPE_DAM (N50 , N50) 

REAL SHAPE_NODE (N50 , N2 0,6) 

REAL XN(N20,N20, 6) , VN (N2 0 , N20 , 6 ) 

J REAL FGEAR (N20,N20, 6) 

— REAL AGEAR (N20 , N20 , N10) , AGSTIF (N20 , N120) 

REAL STIF_GEAR(n20,nl20) 

REAL TMR(N20 ,3,3) 

_ REAL ACUR ( N2 0 ) 

€*****************************************************************0 
INTEGER NGEAR 

— INTEGER INODE, JNODE , ICON, JCON, IROT 

INTEGER AXE, RAD 


REAL PHIN, PS I, GAMA, ALPHA 


REAL SX, SV, FG(6) , FX , FY , FZ , FN , D LENGTH 
REAL F_RADIAL, F_TANG , F AXIAL 
REAL Rl, R2 

REAL GEAR_STIF 

REAL DX, DY, DROT 

REAL XI (6) , VI ( 6 ) , X2(6), V2(6) 

- REAL X21 (6) , V21(6) 

Z* *************************************************************** *c 
:* INITIALIZE VARIABLES *C 

★Tfc***************************************************************^ 

DO 17 IMODE=l , NMODE ( IROT) 

DO 17 ICOORD=l , 6 

_ FGEAR ( IROT , IMODE , ICOORD) = 0.0 

17 CONTINUE 

■ ^*****************************************************************£ 
LC* START CALCULATIONS *C 

C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Q 


NGEAR= CLASS_TABLE (IROT, GEAR) 

DO 100 IGEAR= 1, NGEAR 

LJ DO 10 ICOORD=l,6 

~ XI (ICOORD) = 0.0 

VI (ICOORD) = 0.0 
H- 10 CONTINUE 


C* ******************************************** *c 


H Rl= AGEAR ( IROT , IGEAR , 1 ) 

C* ALL ANGLES IN THE INPUT WERE IN THE DEGREES 
uC* TRANSFORM ALL ANGLES INTO RADIANS 

S DO 11 1=2,5 

AGEAR ( IROT , IGEAR , I) = Pi * AGEAR ( IROT , IGEAR , I ) / 18 0 . 0 
^11 CONTINUE 

~~ PHIN= AGEAR ( IROT , IGEAR , 2 ) 

PSI= AGEAR ( IROT , IGEAR , 3 ) 

|i GAMA= AGEAR ( IROT , IGEAR , 4 ) 

“ ALPHA= AGEAR ( IROT , IGEAR , 5 ) 

£*******************************£ 

gj INODE= NODE_TABLE ( IROT , GEAR, IGEAR) 

DO 20 IMODE=l, NMODE (IROT) 

^ DO 20 ICOORD= 1,6 

s XI (ICOORD) = XI (ICOORD) + 

& XN ( IROT , IMODE, ICOORD) *SHAPE_NODE (INODE, IMODE, ICOORD) 
VI (ICOORD) = VI (ICOORD) + 

a* & VN(IROT, IMODE, ICOORD) *SHAPE NODE ( INODE , IMODE , ICOORD) 
~2 0 CONTINUE 

_ NCON= NODE_CON (INODE, 1) 

LOOP FOR ALL CONNECTIONS ] *C 
— DO 110 ICON= 2, 2*NCON+2 , 3 

* C * FIND CONNECTION 


JNODE= NODE_CON ( INODE , I CON ) 



JCON= NODE_CON ( INODE, ICON+1) 

JTYPE= NODE_CON ( JNODE , ICON+2 ) 

!C ITYPE_ORG = NODE_CON (JNODE , ICON+3 ) 

U ITYPE= JTYPE 

JROT= CLASS_NODE (JNODE) 

:C WRITE (6 6,*) INODE, JNODE, JCON, ITYPE 

C* WORK HERE : CHECK IT AGAIN YOU CAN PUT IT BEFORE BIG LOOP 
NGEAR2= CLASS_TABLE ( JROT , GEAR) 

DO 25 JGEAR=1 , NGEAR2 

— IF ( JNODE .EQ. NODE_TABLE (JROT, GEAR, JGEAR) ) 

& R2= AGEAR(JROT, JGEAR, 1) 

25 CONTINUE 

~C* AS A REFERNCE GEAR IS TAKEN THE ONE WITH SMALLER RADIUS 
ANGLE= A CUR ( I ROT) 

— IF ( R1 .GT. R2 ) ANGLE= ACUR(JROT) 

* :* FIND STIFFNESS FOR THIS CONNECTION 
-C* WORK HERE! 

C* *************************************************************** *q 


GEAR_STIF=0 . 0 
~ NP= NPOINT( ITYPE) 

' C WRITE ( 66 , * ) ' ITYPE= , ' , ITYPE , ' ANGLE= ' , ANGLE 

__C PRINT * , 'BEFORE GSTIF' 

CALL GSTIF ( ITYPE , NP , AGSTIF, ANGLE , STIF_GEAR, GEAR_STIF , 

& ITYPE_DAM, JTYPE) 

C IF (ITYPE .NE. JTYPE) ITYPE - JTYPE 

■^.C PRINT * , 'AFTER GSTIF' 

CC WRITE ( 66 , *) ' GEAR_STIF= ', GEAR_STIF 

"C IF (ITYPE .NE. JTYPE) NODE_CON (JNODE , ICON+2 ) =ITYPE 

******************************** ********************************* 0 . 
X* INITIALIZE VARIABLES *C 

Q* **************************************************************** Q 

• DO 30 ICOORD=l,6 

L X2 (ICOORD) = 0.0 

X21 (ICOORD) = 0.0 
V2 (ICOORD) = 0.0 
V21 (ICOORD) = 0.0 
—3 0 CONTINUE 

*-_* ************** ************************************************* *c 
t .* FIND COORDINATES OF CONNECTED GEAR IN ITS COORD. SYSTEM *C 

JZ * ***************************************** ***********************0 

DO 40 IMODE= 1 , NMODE (JROT) 

L DO 40 ICOORD= 1,6 

X2 (ICOORD) = X2 (ICOORD) + 

& XN ( JROT, IMODE, ICOORD) *SHAPE_NODE (JNODE, IMODE, ICOORD) 

-X V2 (ICOORD) = V2 (ICOORD) + 

& VN (JROT , IMODE , ICOORD) *SHAPE_NODE (JNODE , IMODE , ICOORD) 

40 CONTINUE 

* ************************************************************ ****c 


X:** ************************************** *************************c 
C* COORDINATE TRANSFORMATION OF J-GEAR'S SHAFT COORDINATES INTO *C 
-* I -GEAR'S SHAFT COORDINATES *C 

•^_** ********************************************************* ******C 
C* TRANSFORMATION OF X, Y, Z COORDINATES. ANGLES ARE INVARIANT TO 
n* COORDINATE TRANSFORMATION. 



DO 50 ICOORD= 1,6, 2 
SX= 0.0 
SV= 0.0 

_ DO 51 JCOORD= 1,6, 2 ' 

SX= SX + TMR(JTYPE, ICOORD, JCOORD) *X2 (JCOORD) 
r C SV= SV + TMR(JTYPE, ICOORD, JCOORD) *V2 (JCOORD) 

51 CONTINUE 

X21 (ICOORD) = SX 
CC V21 (ICOORD) = SV 

50 CONTINUE 

* ANGLES REMAIN UNCHANGED, I.E. X21 (angles coord) = X2 (angles coord) 
DO 52 ICOORD= 2,6, 2 ' 

X21 (ICOORD) = X2 (ICOORD) 

52 CONTINUE 

CC WRITE (66 , *) 'DIF ' , X2 ( 1) -X2 1 ( 1) , X2(3)-X21(3) 

****************************************************.j t ^ c ,|| f ,£.j| r .k,kj| f ^^^^^^ 

CALCULATE THE SPRING CONTRACTION ( EXPANSION ) *C 

C****************** ********************************************** * c 

DX= XI (1) -X2 1(1) 

DY= XI ( 3 ) -X2 1(3) 

DROT= R1*X1 (6) -R2*X21 (6) 

_ DLENGTH= ( DX*COS(PHIN) + DY*SIN(PHIN) ) + 

& DROT*COS (PHIN) 

C WRITE (66,*) 'IROT ',IROT,' DL= ' , DLENGTH 

— C WRITE (66, *) 'DX ', DX, ' DY ' , DY, ' DROT ' , DROT 

DVELOSITY= ( Vl(l)-V21(l) )*COS(PHIN) + 
x & ( VI ( 3 ) -V2 1(3) ) *SIN (PHIN) + 

-X & ( Rl*SPEED(IROT) -R2*SPEED(JROT) )*COS(PHIN) 

~***************************ic-k************************************C 

* CALCULATE THE GEAR FORCE AT CONTACT POINT *C 

t** *********************************************************** 

C GEAR_STIF ( JCON) 
w FN= -GEAR_STIF* DLENGTH 

^C WRITE (66 , *) ' IROT= ',IROT,' FN= ' , FN, 

C WRITE (66,*) ' GEAR_ST ' , GEAR_STIF 

C****************************************************************^p 
FX=0 . 0 
’ ' FY=0 . 0 

“ FZ=0 . 0 

DO 12 1=1,6 
12 FG (I ) =0 . 0 

C*************************************************************£££££ 

r* PHIN - PRESSURE ANGLE 
i * PSI - HELIX ANGLE 
w* GAMA - PITCH ANGLE 
C* ALPHA - CONTACT POINT ANGLE 

* F_TANG= FX 
C* F_RADIAL = FY 
C* F_AXIAL= FZ 

CONTACT POINT COORDINATE SYSTEM TRANSFORMATION *C 
F_RADIAL= FN*SIN (PHIN) 

F_TANG= FN*COS (PHIN) *COS (PSI) 

F_AXIAL= FN*SIN (PHIN) *SIN (PSI) 


FX= F_RADIAL 
FY= F_TANG 
_ FZ= F_AXIAL 

0* IN PLANE YZ COORDINATE TRANSFORMATION *C 
AXE= 1 

— RAD= 1 

CALL PLANE ( FX, FY, FZ , GAMA, AXE, RAD ) 

:* IN PLANE XY COORDINATE TRANSFORMATION *C 
AXE= 3 
“ RAD= 1 

CALL PLANE ( FX, FY, FZ , ALPHA, AXE, RAD ) 

_ FG(1)= FX 

FG(2) = FG ( 1) 

FG ( 3 ) = FY 
FG ( 4 ) = FG ( 3 ) 

— FG ( 5 ) = FZ 

^* ********************************************* ****************** *c 

~ DO 70 IM0DE=1 , NMODE ( IROT ) 

DO 70 IC00RD=1 ,5,2 

FGEAR ( IROT , IMODE , ICOORD) = FGEAR(IROT, IMODE, ICOORD) + 

_ & SHAPE_NODE (INODE, IMODE, ICOORD) *FG( ICOORD) 

70 CONTINUE 

DO 75 IMODE=l, NMODE (IROT) 

U FGEAR (IROT, IMODE, 6) = FGEAR (IROT, IMODE, 6) + 

& SHAPE_NODE (INODE, IMODE, 6) *FY*R1 
75 CONTINUE 

— 110 CONTINUE 
100 CONTINUE 

JC WRITE (66 , *) ' FGEAR= ' , FGEAR ( IROT, 1, 1) 

RETURN 
m END 

************************************************************0 

c***********************************************************c 

jC * ************************************** Q 

| J* CODE FOR GEAR STIFFNESS CALCULATION *C 
»^***************************************£ 

SUBROUTINE GSTIF( ITYPE, NPOINT, AGSTIF, ANGLE, 

& STIF_GEAR, GEAR_STIF, ITYPE_DAM, JTYPE ) 

“ PARAMETER N10=10 

PARAMETER N20=20 
*4 PARAMETER N50=50 

g PARAMETER N12 0=120 

INTEGER ITYPE, NPOINT, JTYPE 
INTEGER ITYPE_DAM (N50 , N50) 

REAL AGSTIF (N20,N120) , ANGLE 
REAL STIF_GEAR(n20,nl20) 

REAL GEAR_STIF 
■c***********************c 
INTEGER I POINT 

_ REAL A, B, ARATIO, SDIF, AMAX 

'!** F*** REWRITE WITH FORTRAN FUNCTIONS : MOD & ETC. *C 
C WRITE (66 , *) ' ANGLE= ANGLE 



ITEETH = 1 

10 IF ( ANGLE .GE. 360.0 ) THEN 

ANGLE= ANGLE - 360.0 
ITEETH = 1 
GOTO 10 
ENDIF 

20 IF ( ANGLE .GE. 0.0 ) GOTO 25 
ANGLE= ANGLE + 360.0 
ITEETH = 1 
GOTO 20 


25 CONTINUE 

AMAX= AGSTIF ( ITYPE , 1 ) 

30 IF ( ANGLE .LE. AMAX ) GOTO 40 
ANGLE= ANGLE - AMAX 
ITEETH - ITEETH +1 
GOTO 30 


40 CONTINUE 

r C IF ( ITEETH .GE. 28) PRINT * , 7 ITEETH & ANGLE= 

IF ( ITYPE_DAM(JTYPE, ITEETH) .NE. ITYPE) THEN 
ITYPE = ITYPE DAM (JTYPE, ITEETH) 

END IF 


ITEETH, ANGLE 


^,it it it it it it it it it it it it ik it it it Q 

if (iteeth .EQ. 10) then 

■c write (6 6,*) 'itype 7 , itype, 'iteeth 7 , iteeth, 7 j type' , jtype 

C write (66,*) 'angle 7 , angle J 

:-a write (66, *) 'stif gear 7 , (stif gear (itype, i),i=l,100) 

2 end if - ' 

GEAR_STIF= STIF_GEAR ( ITYPE , 1 ) 


; * THERE IS A BUG = FIND & FIX= RUN 

— DO 50 IPOINT= 2 , NPOINT 

IF ( ( ANGLE .GE. AGSTIF (ITYPE, 3+IPOINT) ) 

- - & . AND « 

& ( ANGLE .LT. AGSTIF (ITYPE, 4+IPOINT) ) ) THEN 


A= ANGLE - AGSTIF (ITYPE, 3+IPOINT) 

B= AGSTIF(ITYPE, 4+IPOINT) - AGSTIF (ITYPE, 3+IPOINT) 
ARATIO= A/B 

SDIF= STIF_GEAR( ITYPE, I POINT) - STIF GEAR (ITYPE, IPOINT- 
GEAR_STIF= STIF_GEAR ( ITYPE, IPOINT-1) + SDIF*ARATIO 


1 ) 


ENDIF 


-50 CONTINUE 

w 

WRITE (66,*) 7 ITEETH= 7 ITEETH 

ITYPE = JTYPE 
w RETURN 

END 


f* *************************************** ************************ 

C* *************************************************************** 

* SUBROUTINE FOR IN-PLANE COORDINATE TRANSFORMATION 
«»»* ************************************* * * * * * * * * * * * * * * * * * * * * * * * * * * 


rm 

|=3 


SUBROUTINE PLANE ( X, Y, Z, ANGLE, AXE, 

PARAMETER PI=3 . 14 

INTEGER AXE, DIR, RAD 

REAL X, Y, Z, ANGLE 

REAL OLD_A , OLD_B, OLD C 

REAL NEW_A, NEW_B, NEW_C 


RAD ) 


*C 

*C 

*c 

*c 


IF ( RAD .EQ. 0 ) ANGLE = PI*ANGLE/180 . 0 


DIR= IABS (AXE) /AXE 


IF ( AXE .EQ. 1 ) THEN 
OLD_A= Y 
_ OLD_B= Z 

OLD_C= X 
ENDIF 

— IF ( AXE .EQ. 2 ) THEN 

OLD_A= X 

OLD_B= Z 
OLD_C= Y 

— ENDIF 

IF ( AXE .EQ. 3 ) THEN 

OLD_A= X 
” OLD_B= Y 

OLD_C= Z 
ENDIF 

NEW_A= OLD_A* COS (ANGLE) - OLD_B* SIN (ANGLE) 

NEW_B= OLD_A*SIN (ANGLE) + OLD_B*COS (ANGLE) 

NEW_C= DIR*OLD_C 

IF ( AXE .EQ. 1 ) THEN 
Y= NEW_A 
Z= NEW_B 
~ X= NEW_C 

ENDIF 

IF ( AXE .EQ. 2 ) THEN 
X= NEW_A 
Z= NEW_B 
Y= NEW_C 

— ENDIF 

IF ( AXE .EQ. 3 ) THEN 
X= NEW_A 
“ Y= NEW_B 

Z= NEW_C 
ENDIF 

“ RETURN 

END 

i’*******************************************************************C 



"C * *************************************************************** ******* c 
C************************ •• THE FINAL INSULT " *************************£ 
3*** ********************************************************** 

SUBROUTINE LAT (NUNIT, NCOUNT, NROT, NSTAT , NBEAP, LBEAR, STIF2 
& DAMP2,STIFA, FREQ, SHAPE, NMODE,NCON BEAR, NSWITCH) ' 

PARAMETER NROT_SIZE=10 “ 

PARAMETER NSTAT_SIZE=100 
PARAMETER N100=100 
PARAMETER N50=50 
PARAMETER N20=20 
PARAMETER N10=10 
PARAMETER N5=5 


INTEGER NMODE (N50) 

INTEGER NM (N50 ) 

INTEGER NROT 

INTEGER NUNIT (N20) , NSTAT (NROT_SIZE) 

INTEGER NBEAR(NROT SIZE) , LBEAR (NROT SIZE, 10) 
INTEGER NSWITCH (N10) ~ 


REAL SHAPE (N50,N50,N50,4) 

REAL FREQ (N10 , N10 , NIO) 

REAL LMODE (NROT_SIZE , N10 , NSTAT SIZE) 
REAL SLOPE (NROT_SIZE, N10 , NSTAT_SIZE) 


INTEGER N, NB, NNCT 
INTEGER NBC 
REAL SPI,SPL, DSP, DDIN 
REAL CRT(NIO) 

REAL TCRT (N10 , N5) 

REAL DDPC (N10 , N100) , EEYTH (N10 , N100) 

REAL TDMY (N10 , N10 , N5) , TEEYTH (N10 , N100 , N5) 

REAL TW(N100 , 5) , TWMOD (N50 , N5) 

REAL AKK(NIO) ,TAKK(N10,N5) , AKRR(NIO) ,TAKR(N10,N5) 

REAL DEFL(NIOO) , EYTH(NIOO), WMOD(NIOO), TMX(N100,N100) 
REAL EYI(NIOO) , EY2 (N100) 

REAL DPC(NIOO) , EANI(NIOO), EAN2(N100) 

REAL KXX(NIO) , KYY (N10) , KRX (N10) , KRY(NIO) 

REAL TDDPC (N10,N100,N5) 

REAL SHAFT_WEIGHT(NSTAT_SIZE) ,TOTAL_WEIGHT, W (NSTAT SIZE) 
REAL EXTERNAL__WEIGHT (NSTAT SIZE), WEIGHT STAT (NSTAT SIZE) 
REAL DOUT (NSTAT_SIZE) , DIN (NSTAT_SIZE) - 

REAL ITRANS (NSTAT_SIZE) , IPOLAR (NSTAT SIZE) 

REAL DX(NSTAT_SIZE) ,RO(NSTAT_SIZE) 

REAL El (NSTAT_SIZE) , INERTIA (NSTAT SIZE) 

REAL EM (NSTAT_SIZE) - 

REAL STIF_BEAR (NROT_SIZE , N10 ,6,6) 

REAL DAMP_BEAR (NROT_SIZE , N10 , 6 , 6) 

REAL STIF2 (NROT_SIZE , N10 , 6 , 6 ) , STIFA(NROT SIZE,N10, 6, 6) 
REAL DAMP2 (NROT_SIZE , N10 , 6,6) ~ 


******************************************************** *c 
REAL WJ(NSTAT_SIZE) ,GJ (NSTAT SIZE) 

REAL GGO(NSTAT_SIZE) ,TK (NSTAT SIZE) 

: REAL AT (N5 , N5) ” 


REAL FFR(NIO) , FV(N100 ,N10) , TV (N100 , N10) , WF (N100 , N100) 

REAL FR (N10 , N5) , FFV (N10 , N100 , N5) 

REAL ZF(N10,N5) , ZFV(N10 , N100 , N5) 

REAL SP1, SP2 , SP3 

s * ************************* ***************************************c 


C** FOR SYSTEM NATURAL FREQ. CALCULATIONS ******************C 

= REAL AK , AKR , STIF_AV (N10 , N10 , 4 , 4 ) 

M REAL FMODE(N10,N10,4) , SBARK, SB ARC 


1 



REAL KBAR (N10 , N10) , CBAR (NIO , NIO) 

REAL K2BAR(N10,N10) , C2BAR (NIO , NIO) 

REAL AGLOB ( 10 ,24,24) , ASUPER (N100 , N100 ) 

— INTEGER NCON_BEAR (NIO , NIO, NIO, NIO) 

C************************* INPUT DATA *************************** c 

^* ****************************************************************£ 
l************************* START MAIN BODY ***********************c 

IR0T=1 

READ (NUNIT (1) , *) , 

_ READ (NUNIT (1),*) NROT 

C********** START ENTERING SHAFT AND BEARINGS DATA ***************c 
10 CONTINUE 

READ (NUNIT (1),*) 

READ (NUNIT (1),*) 

READ (NUNIT (1) , * ) NSTAT(IROT) 

READ (NUNIT (1),*) 

— READ (NUNIT (1) , *) ( EXTERNAL_WEIGHT (J) , DX (J) , DOUT (J) , DIN (J) , 

+ IPOLAR(J) ,ITRANS(J) ,EM(J) ,RO(J) ,GGO(J) ,TK(J) , J=1 , NSTAT (IROT) ) 
; READ (NUNIT (1),*) 

READ (NUNIT ( 1) , *) 

— READ (NUNIT (1) , * ) NBEAR(IROT) 

READ (NUNIT (1),*) 

READ (NUNIT (1) , *) ( LBEAR ( IROT, IBEAR) , IBEAR=1 , NBEAR ( IROT) ) 

• • DO 20 IBEAR=1 , NBEAR ( IROT) 

READ (NUNIT (1),*) 

READ (NUNIT (1),*) 

READ (NUNIT (1) ,*) (( STIF BEAR (IROT, IBEAR, I , J) , J=1 , 6) , 1=1 , 6) 

L_ READ (NUNIT (1),*) 

READ (NUNIT (1) , *) ( ( DAMP BEAR(IROT, IBEAR, I, J) , J=l, 6) , 1=1, 6) 

■ 20 CONTINUE “ 

READ (NUNIT (1) , *) 

— READ (NUNIT ( 1 ) , * ) 

READ (NUNIT (1) , *) NMODE ( IROT) , SPI , SPL, DSP 
READ (NUNIT (1),*) 

READ (NUNIT (1) , *) NBC 

CCC NMODE ( IROT ) =3 

^***** *************** ************* ******* * **************** * ***c 

C*** INITIALIZE MODE SHAPES AND FREQUENCIES *****************0 
_ DO 70 ISPEED=1, NMODE (IROT) 

FREQ (IROT, ISPEED, 1) = 0.0 

— DO 70 ISTAT=1, NSTAT (IROT) 

SHAPE (IROT, ISPEED, ISTAT, 1)= 0.0 

M SHAPE (IROT, ISPEED, ISTAT, 2) = 0.0 

^70 CONTINUE 


c*************************************************************c 

* ************************** 

S* * CALCULATE EFFECTIVE INERTIA MOMENT OF ROTORS * 

C* ************************** 




So 


PI=3. 14159 
E=EM (1) 

DO 30 1=1, NSTAT (IROT) 

INERTIA (I) = PI*( DOUT ( I ) * * 4 - DIN(I)**4 ) /64 
El (I) = E* INERTIA (I) 

SHAFT_WEIGHT ( I ) = PI*( DOUT(I)**2 - DIN(I)**2 ) *DX (I) *RO (I) /4 


W(l)= SHAFT_WEIGHT(l)/2.0 + EXTERNAL WEIGHT (1) 
TOTAL_WEIGHT= W(l) 

TOTAL_LENGTH= DX(1) 


DO 40 1=2, NSTAT (IROT) 



WEIGHT_STAT ( I ) = SHAFT_WEIGHT ( 1-1) /2 . 0 + SHAFT WEIGHT (I) /2 . 0 
& + EXTERNAL WE IGHT ( I ) 

W(I)= WEIGHT_STAT ( I ) , 

- TOTAL_WEIGHT= TOTAL_WE IGHT + WEIGHT_STAT (I) 

TOTAL_LENGTH= TOTAL LENGTH + DX(I) 

40 CONTINUE 

- I POLAR ( 1) = IPOLAR(l) + INERTIA (1) *RO(l) *DX(1) 

ITRANS (1) = ITRANS ( 1) + 

& SHAFT WEIGHT (1) * ( ( DOUT(1)**2.0 + DIN(1)**2.0 )/16.0 + 

_ & (( DXX 1 ) / 2 . 0 ) **2 . 0 )/3.0 )/2.0 

DO 50 1=2 , NSTAT ( IROT) 

IPOLAR(I) = I POLAR (I) + RO ( I ) *INERTIA ( I ) *DX (I) + 

- & INERTIA(I-l) *DX(I-1) *RO(I-l) 

ITRANS (I) = ITRANS (I) + 

& SHAFT WEIGHT (I) * ( (DOUT (I) **2 . 0+DIN (I)**2.0)/16.0 + 

& (( DXTl)/2.0 ) **2 . 0 )/3.0 )/2 . 0 + 

- & SHAFT_WEIGHT (1-1) * ( (DOUT (1-1) **2 . 0+DIN (I — 1 ) **2 . 0) /16 . 0 + 

& (( DX(I-l)/2.0 ) **2 . 0 )/3.0 )/2.0 

50 CONTINUE 

"£*************** OUTPUT SHAFT DATA CALCULATIONS ********************0 
WRITE (NUNIT (NCOUNT) ,*) 

WRITE (NUNIT (NCOUNT) , *) 7 DATA FOR ROTOR 7 # 7 , IROT 

_ WRITE (NUNIT (NCOUNT) , *) 

WRITE (NUNIT (NCOUNT) , *) 7 I W(I) DX(I) DOUT(I) DIN(I) INERTIA ( I ) 

1 I POLAR (I) ITRANS (I) El (I) 7 

WRITE (NUNIT (NCOUNT) , *) 

- WRITE (NUNIT (NCOUNT) ,61) ( I , W (I) , DX ( I) , DOUT (I) , DIN (I ) , INERTIA ( I) , 

1 IPOLAR(I) ,ITRANS(I) ,EI(I) , 1=1 , NSTAT ( IROT) ) 

WRITE (NUNIT (NCOUNT) ,*) 

WRITE (NUNIT (NCOUNT) ,*) 'TOTAL WEIGHT 7 , 'TOTAL LENGTH 7 
~ WRITE (NUNIT (NCOUNT) , *) TOTAL_WE I GHT , TOTAL LENGTH 

WRITE (NUNIT (NCOUNT) , *) 

DO 621 IBEAR=1 , NBEAR ( IROT) 

_ KXX (I BEAR) = STIF_BEAR(IROT, I BEAR, 1,1) 

KYY (I BEAR) = STIF_BEAR ( IROT, I BEAR, 3,3) 

KRX (IBEAR) = STIF_BEAR ( IROT , IBEAR ,2,2) 

KRY ( IBEAR) = STIF_BEAR ( IROT , IBEAR ,4,4) 

- WRITE (NUNIT (NCOUNT) ,*) 'BEARING STIFFNESS IN X- AND Y- DIRECTIONS' 

WRITE (NUNIT (NCOUNT) ,*) 7 KXX ', 7 KYY 7 

WRITE (NUNIT (NCOUNT) , *) KXX (IBEAR) , 7 7 ,KYY (IBEAR) 
y WRITE (NUNIT (NCOUNT) , *) 7 AVERAGE SPRING BEARING STIFFNESS 7 

WRITE (NUNIT (NCOUNT) , *) 0 . 5* (KXX (IBEAR) + KYY (IBEAR) ) 

WRITE (NUNIT (NCOUNT) , *) 7 AVERAGE ROTATIONAL BEARING STIFFNESS 7 

WRITE (NUNIT (NCOUNT) , *) 0 . 5* (KRX (IBEAR) + KRY (IBEAR) ) 

_621 CONTINUE 

WRITE (NUNIT (NCOUNT) , *) 7 DATA FOR LATERAL MODE CALCULATIONS: 7 

SPI=INITIAL SPEED, SPL=FINAL SPEED, DSP=SPEED INCREMENT-RPM *C 
= WRITE (NUNIT (NCOUNT) , *) 7 SPI 7 , 7 SPL DSP 7 

WRITE (NUNIT (NCOUNT) , *) SPI, SPL, DSP 

WRITE (NUNIT (NCOUNT) , *) 7 TYPE OF BOUNDARY CONDITIONS 7 

WRITE (NUNIT (NCOUNT) , *) NBC 

C********************* PART TWO ******************************0 
**************************************************************£ 


*********************** 

C * CALCULATE ROTOR'S LATERAL MODAL SHAPE * 

r *********************** 


C* NC=LOCAL CRITICAL SPEED 'NO. *C 
NC=1 

MA=0 ! 



MB=0 

LN=1 

LN=LN+3 

- DDIN=DSP 
SPD=SPI 
DETP=0 . 

X******** START LOOP 290-610 *********C 
290 1=1 

J=1 

_ SPSQ=SPD*SPD 

ANSP=SPD*0. 10471976 
ANSP2=ANSP*ANSP 1 

VP=0. 1 

_ ZMP=0 . 

EYP=0 . 

ETHP=1 . 0 
M=1 

^00 1=1+1 

11=1-1 

IF (II-LBEAR(IROT, J) ) 330,310,330 
10 AK= (KXX ( J) +KYY ( J) ) /2 . 0 

AKR= (KRX ( J) +KRY ( J) ) /2 . 

AKK(J) =AK 
AKRR(J) =AKR 

TAKK(J, IROT) =AKK ( J) 

TAKR ( J , IROT ) =AKRR ( J ) 

IF (J-NBEAR(IROT) ) 320,340,340 
20 J=J+1 

- GO TO 340 

330 AK=0 . 0 

AKR=0 . 0 

40 VP=VP+(W(I-1) *ANSP2/386 . 4-AK) *EYP 

- ZMP=ZMP+AKR*ETHP-ANSP2* (ITRANS (1-1) )*ETHP/386.4 
EY=E YP+DX ( I - 1 ) *ETHP+DX (I — 1 ) **2*ZMP/ (2 . E6*EI (1-1) ) + 

1 DX(I-l) **3*VP/(6.E6*EI(I-1) ) 

U ETH=ETHP+DX ( I- 1 ) *ZMP/ ( 1 . E6*EI (1-1) ) + 

1 DX(I-l) **2*VP/(2.E6*EI(I-1) ) 

ZM=ZMP+DX ( 1-1) *VP 
I V=VP 

U IF (M.EQ.2) GO TO 350 

EY1 (I) =EY 
EAN1 (I) =ETH 

IF (I.GT.NSTAT(IROT) ) GO TO 360 
“ ZMP=ZM 

VP=V 
EYP=EY 

_ ETHP=ETH 

GO TO 300 
150 EY2 (I) =EY 

EAN2 ( I ) =ETH 

- ZMP=ZM 
VP=V 
EYP=EY 
ETHP=ETH 

~ IF (I . GT. NSTAT ( IROT) ) GO TO 370 

GO TO 300 
60 M=2 

_ ZM1=ZM 

VR1=V 
J=1 
1=1 

- EYP=1 . 

ZMP=0 . 

ETHP=0 . 

VP=0. 



“ GO TO 300 

370 DET=VR1*ZM-V*ZM1 

IF (ABS(DETP) .LT. 0.0001) GO TO 420 
_ IF (MA.EQ.l) GO TO 400 

IF (ABS(DET) .LT. 1. ) GO TO 450 
IF (DETP*DET) 380,420,420 
180 DOLD=DETP 
— 190 MA=1 

IF (ABS(DET) .LT. 1. ) GO TO 450 
IF (DDIN.LT. l.E-6) GO TO 450 
DDIN=DDIN/2 . 

— DETPP=DETP 
DETP=DET 
SPD=SPD-DDIN 
GO TO 290 

"400 IF (ABS(DET) .LT. 1. ) GO TO 450 
IF (DOLD*DET) 390,420,410 
i =10 CONTINUE 

Id IF (ABS(DET) .LT. 1. ) GO TO 450 

IF (DDIN.LT. l.E-6) GO TO 450 
DDIN=DDIN/2 . 

SPD=SPD+DDIN 

— DETPP=DETP 
DETP=DET 
GO TO 290 

20 IF (LN-54) 440,440,430 

“4 30 CONTINUE 

LN=1 

40 CONTINUE 

LN=LN+1 
SPD=SPD+DSP 
DDIN=DSP 

IF(NC.GT.NMODE(IROT) ) GO TO 610 

— IF (SPD.GT.SPL) GO TO 610 
DETPP=DETP 

. DETP=DET 

p SSPD=SPD 

**• GO TO 290 

450 MA=0 

id LN=LN+1 

& IF (LN-50) 470,470,460 

*460 CONTINUE 

LN=1 

1 70 CONTINUE 

w CRT (NC) =SPD 

TCRT ( NC , I ROT ) =CRT ( NC ) 

^** ******************************** 

» NC=NC+1 

* LN=LN+3 

EY1 (1) =0 . 

|~j EY2 (1) =1 . 

= DTX=0 . 

1=1 

IF (LN-50) 490,490,480 
80 CONTINUE 

— LN=1 

490 CONTINUE 

LN=LN+2 

00 DEFL(I) =V*EY1 (I) -VR1*EY2 (I) 

— IF (I.NE.l) GO TO 510 
EYTH ( I ) =V 

GO TO 520 

10 EYTH(I) =EAN1 (I) *V-EAN2 (I) *VR1 
1>2 0 DEFA=ABS ( DEFL ( I ) ) 

DMXA=ABS (DTX) 

— 1 = 1+1 



~ IF (DEFA-DMXA) 540,540,530 

530 DTX=DEFL(I— 1) 

>40 IF ( I-NSTAT ( IROT) ) 550,550,560 
_>50 GO TO 500 

560 DO 570 1=1, NSTAT( IROT) 

DPC ( I ) =DEFL ( I ) /DTX 
EYTH ( I ) =EYTH ( I ) /DTX 
_ EEYTH (NC- 1,1) =EYTH ( I ) 

570 DDPC (NC-1 , I ) =DPC ( I ) 

DO 600 1=1 , NSTAT ( IROT) 

LN=LN+1 

— IF (LN-54) 590,590,580 

580 CONTINUE 

LN=1 

>90 CONTINUE 

LN=LN+1 

600 CONTINUE 

SPD=SSPD+DSP 
_ DETP=0 . 

GO TO 290 
610 CONTINUE 

;************** END OF THE LOOP 290-610 *************************C 
-j ** ************* ******** ************************************** ***c 
DO 6200 IBEAR= 1 , NBEAR ( IROT) 

DO 6201 1=1,6 
DO 6201 J=1 , 6 

STIF2 ( IROT , IBEAR , I , J ) = STIF_BEAR(IROT, IBEAR, I, J) 

DAMP2 ( IROT , IBEAR , I , J) = DAMP_BEAR ( IROT , IBEAR , I , J) 

6201 CONTINUE 
^ 6200 CONTINUE 

c****************************************************************c 
C* SUBSTRACTING AVERAGE BEARING STIFFNESS FROM STIFFNESS MATRIX *C 
DO 620 I BEAR=1, NBEAR (IROT) 

_ STIF_AV ( IROT , IBEAR ,1,1) =0 . 5 * ( STIF_BEAR (IROT, IBEAR, 1 , 1) + 

& STIF_BEAR(IROT, IBEAR, 3,3) ) 

STIF_AV (IROT, IBEAR, 2,2) =0 • 5* ( STIF_BEAR(IROT, IBEAR, 2 , 2) + 

& STIF_BEAR(IROT, IBEAR, 4,4) ) 

— AK= STIF_AV (IROT, IBEAR, 1, 1) 

AKR= STIF AV( IROT, IBEAR, 2, 2) 

STIF_BEARllROT, IBEAR, 1, 1) =STIF_BEAR(IROT, IBEAR, 1, 1) -AK 
_ STIF_BEAR ( IROT , IBEAR ,3,3) =STIF_BEAR ( IROT , IBEAR, 3,3) -AK 

STIF_BEAR ( IROT , IBEAR ,2,2) =STIF_BEAR ( IROT , IBEAR, 2,2) -AKR 
STIF_BEAR ( IROT , IBEAR, 4,4) =STIF_BEAR ( IROT , IBEAR, 4,4) -AKR 
“20 CONTINUE 

**************************************** ********************** *c 
C* FOR THE TRANSIENT PART OF THE PROGRAM *C 

ri DO 6210 IBEAR= 1 , NBEAR ( IROT) 

— DO 6211 1=1,4 
DO 6211 J=1 , 4 

“ STIFA ( IROT , IBEAR, I , J) = STIF_AV (IROT, IBEAR, I, J) 

6211 CONTINUE 
W 6210 CONTINUE 

DO 6220 IBEAR= 1 , NBEAR ( IROT) 

— DO 6221 1=5,6 
DO 6221 J=5 , 6 

= STIFA (IROT, IBEAR, I , J) = 0.0 

^ 6221 CONTINUE 
^6220 CONTINUE 

c*************************************+**************************c 

AK=0 . 0 

= AKR=0 . 0 

DO 60 I BEAR=1, NBEAR (IROT) 

WRITE (NUNIT (NCOUNT) , *) ' STIFFNESS MATRIX FOR ' , IBEAR, ' BEARING ' 

WRITE (NUNIT (NCOUNT) ,2) ( ( STIF2 (IROT, IBEAR, I, J) , J=l, 6) , 1=1, 6) 



WRITE (NUNIT (NCOUNT) , *) 

WRITE (NUNIT (NCOUNT) , *) ' DAMPING MATRIX FOR ' , IBEAR , ' BEARING ' 

WRITE (NUNIT (NCOUNT) ,2) ( ( DAMP2 (IROT, IBEAR, I, J) , J=l, 6) , 1=1, 6) 
U 2 FORMAT (2X, 6F12.3 ) 

60 CONTINUE 


;*********************************************************** * c 
-C** NCS -> NUMBER OF CRITICAL SPEED= NSPEED II=ISPEED 
NSPEED= NC-1 


!** OUTPUT THE RESULTS OF TIJE ABOVE CALCULATIONS ******************c 

— DO 650 ISPEED=1, NSPEED 
WMOD ( ISPEED) =0 . 

DO 630 ISTAT=1 , NSTAT ( IROT) 

__ 30 WMOD ( ISPEED) =WMOD( ISPEED) + 

& ITRANS (ISTAT) *EEYTH( ISPEED, ISTAT) **2.0 + 

& W ( ISTAT) *DDPC( ISPEED, ISTAT) **2.0 

WMOD ( ISPEED) =WMOD( ISPEED) /38 6. 4 

— TWMOD (ISPEED, IROT) =WMOD( IS PEED) 


DO 640 ISTAT=1, NSTAT (IROT) 

** EEYTH - SLOPE OF THE BEAM 

w EEYTH (ISPEED, ISTAT) =EEYTH( ISPEED, ISTAT) /(WMOD (ISPEED) **0.5) 

C** DDPC - DEFLECTION OF THE BEAM 

DDPC ( ISPEED , ISTAT ) =DDPC ( ISPEED , ISTAT) / (WMOD ( ISPEED) * *0 . 5 ) 

_** REAL LMODE(NROT,NMODE, NSTAT) - LATERAL MODE SHAPES OF THE SYSTEM 
LMODE (IROT, ISPEED, ISTAT) = DDPC (ISPEED, ISTAT) 
r** TEEYTH - SLOPE OF THE BfeAM 

SLOPE (IROT, ISPEED, ISTAT) = EEYTH (ISPEED, ISTAT) 

—** SHAPE ARRAY KEEPS THE RESULTS OF CALCULATIONS **************c 
SHAPE ( IROT , ISPEED , ISTAT , 1 ) = LMODE ( IROT , ISPEED , ISTAT) 

SHAPE(IROT, ISPEED, ISTAT, 2) = SLOPE ( IROT , ISPEED, ISTAT) 

—640 CONTINUE 
650 CONTINUE 


OUTPUT DATA ' 

LATERAL FREQUENCIES AND MODE 


ROTOR 


CC 


/ / ^ / IROT 
' LATERAL FREQUENCIES 


'NO. 


FREQUENCY ( HZ ) ' , 


£>** 

u 


fc_i 

EJ 


WRITE (NUNIT (NCOUNT-2 ), *) ' 

WRITE (NUNIT (NCOUNT-2 ), *) ' 

& SHAPES ' 

WRITE (NUNIT (NCOUNT-2 ), *) 

WRITE (NUNIT (NCOUNT-2 ), *) ' 

WRITE (NUNIT (NCOUNT-2) , *) 

DO 660 ISPEED=1, NSPEED 
WRITE (NUNIT (NCOUNT-2 ), *) 

WRITE (NUNIT (NCOUNT-2 ), *) 

1 ' MODAL WEIGHT 7 

WRITE (NUNIT (NCOUNT-2 ) ,*) ISPEED, CRT ( ISPEED) /60 . 0 , WMOD(ISPEED) 

************************************************************C 

FREQ (IROT, ISPEED, 1 ) =CRT ( ISPEED) 

FREQ (IROT, ISPEED,^ ) =CRT ( ISPEED) 

WRITE (NUNIT (NCOUNT-2 j , *) 

WRITE (NUNIT (NCOUNT-2 ), *) 

WRITE (NUNIT (NCOUNT-2 ), *) 

WRITE (NUNIT (NCOUNT-2 ),*) r #STATION 
WRITE (NUNIT (NCOUNT-2 ) , * ) 

DO 660 ISTAT=1, NSTAT (IROT) 

WRITE (NUNIT (NCOUNT-2) , *) ISTAT, LMODE (IROT, ISPEED, ISTAT) , 

& SLOPE (IROT, ISPEED, ISTAT) 

CONTINUE 


MODE SHAPES FOR ROTOR NO . ' , IROT 

DEFLECTION SLOPE 


WRITE (NUNIT (NCOUNT-2 ) , * ) 

WRITE (NUNIT (NCOUNT-2 ), *) ' CHECK-OUT THE ORTHOGONALITY OF 

& THE MODE SHAPES ' 

WRITE (NUNIT (NCOUNT-2 ), *) ' MODE SHAPES ARE NORMALIZED WITH 

& RESPECT TO MASS MATRIX 7 



WRITE (NUNIT (NCOUNT-2 ) , * ) 

DO 680 ISPEED=1 , NSPEED 
DO 670 JSPEED=1, NSPEED 
~ TMX (ISPEED, JSPEED) =0. ' 

DO 670 ISTAT=1 , NSTAT (IROT) 

TMX ( ISPEED , JSPEED) =TMX ( ISPEED, JSPEED) + 

& W(ISTAT) *LMODE (IROT, ISPEED, ISTAT) * 

— & LMODE (IROT, JSPEED, ISTAT) /386 . 4+ 

& ITRANS (ISTAT) *SLOPE (IROT, ISPEED, ISTAT) * 

& SLOPE (IROT, JSPEED, ISTAT)/386. 4 

_I70 CONTINUE 

WRITE (NUNIT (NCOUNT-2 ), *) ( TMX (ISPEED, JSPEED) ,JSPEED=1, NSPEED ) 
680 CONTINUE 

NCT=NSPEED 

CCCCCCCCC NMODE (IROT) =NCT 

NM ( IROT) =NCT 

— WRITE (NUNIT (NCOUNT) ,*) 

WRITE (NUNIT (NCOUNT) , *) 'NUMBER OF FOUND NATURAL FREQUENCIES' 
WRITE (NUNIT (NCOUNT) , *) NM(IROT) 

IF ( NMODE (IROT) - NM(IROT) .NE. 0 ) THEN 

— WRITE (NUNIT (NCOUNT) ,*) 'THE HIGHER FREQUENCIES AND MODE SHAPES' 

WRITE (NUNIT (NCOUNT) , *) ' WILL BE TAKEN AS ZEROS ' 

ENDIF 

C******************* TORSIONAL PART *******************************0 

** N= NSTAT 

^** NM -> NMODE=NSPEED 

C *********************************************************** 

~ * THIS PROGRAM FOR CALCULATE Z -DIRECTION EIGEN PROBLEM * 

^ *********************************************************** 

C* NTYPE=1 CORRESPONDS TO TORSIONAL VIBRATION 
* NTYPE=2 CORRESPONDS TO Z-DIRECTION ( THRUST ) VIBRATION 

w NTYPE=1 

700 CONTINUE 

! >* NMODE=l 

_** IROT=NROT 

C** NE=IROT i 

= NSTATX=NSTAT ( IROT) ; 

m IF (NTYPE.EQ.l) THEN , 

w DO 710 ISTAT=1, NSTAT (IROT) 

GGO (ISTAT) =GGO (ISTAT) *10000000 . 0 
z.i GJ ( ISTAT ) =3 2 . 0 * DX ( ISTAT ) / ( PI *GGO ( ISTAT) * 

y & ( DOUT (ISTAT) **4-DIN ( ISTAT) **4 ) ) 

WJ ( ISTAT) =PI*RO (ISTAT) *DX(ISTAT) * 

& ( DOUT (ISTAT) **4-DIN ( ISTAT) **4)/32.0/386.4+ 

“ & ITRANS ( ISTAT) /3 8 6. 4 

LJ10 CONTINUE 

ENDIF 

ti IF (NTYPE.EQ. 2) THEN 

DO 720 ISTAT=1, NSTAT (IROT) 

65 TK( ISTAT) =0.0 

EM ( ISTAT) =EM( ISTAT) *1000000.0 
“ GJ (ISTAT) = 4 . 0*DX ( ISTAT) / (PI*EM (ISTAT) * 

y & (DOUT (ISTAT) **2-DIN (ISTAT) **2 ) ) 

WJ (ISTAT) = PI*RO (ISTAT) *DX( ISTAT) * 

& ( DOUT ( ISTAT ) **2-DIN (ISTAT) **2 ) /4 . 0/386 . 4+ 

& EXTERNAL WEIGHT ( ISTAT) /3 86 . 4 

i— 72 0 CONTINUE 

ENDIF 

***** ** *************** *** * * * Q 



730 


MCY=0 
ST=20 . 0 
AO 1=0 . 0 
A01=A01+ST 
A02=A01+ST 
c* *********************************** *c 

CALL SBC (NSTATX, NBC , AOl, B1,GJ,WJ,TK) 

CALL SBC (NSTATX, NBC , A02 , B2 , GJ , WJ , TK) 

0************************i****i*i**i**C 

IF ( (B1*B2 ) . LE .0.0) GOTO 740 
GO TO 730 
_ 740 MCY=MCY+1 

IF (MCY . GT . NMODE ( IROT) ) GOTO 770 

A1=A01 

A2=A02 

-750 IF (ABS (A1-A2 ) .LT.0.1) GOTO 760 
AAl=Al+0. 618* (A2-A1) 

]************************^r***************C 
CALL SBC (NSTATX , NBC, AA1 , BB , GJ , WJ, TK) 

"C* ***************************** **********0 
IF ( (B1*BB) .LE. 0.0 ) THEN 
A2=AA1 
„ B2=BB 

ELSE 
A1=AA1 
B1=BB 

— ENDIF 
GOTO 750 

760 FR(MCY, IROT) =0.5* (A1+A2) 

IF ( NTYPE .GT. 1 ) ZF (MCY, IROT) =FR (MCY, IROT) 

- FFR (MCY) =FR (MCY, IROT) /2.0/PI 
GOTO 730 

770 CONTINUE 

N1=NSTAT (IROT) +1 

~ DO 780 IMODE=l, NMODE (IROT) 

FV ( 1 , IMODE ) =1 . 0 

TV ( 1 , IMODE ) =-WJ ( 1 ) *FR( IMODE, IROT) **2 
_ DO 780 ISTAT=2 , NSTAT ( IROT) 

TV (ISTAT, IMODE) =TV ( ISTAT-1 , IMODE) * ( 1.0+ 

& ( TK( ISTAT) -WJ (ISTAT) *FR (IMODE, IROT) **2) *GJ (ISTAT) ) 

& +FV(ISTAT-1, IMODE) * (TK(ISTAT) -WJ (ISTAT) *FR(IMODE, IROT) **2 ) 
~ FV (ISTAT, IMODE) =TV (ISTAT-1 , IMODE) *GJ (ISTAT) + FV ( ISTAT-1 , IMODE) 

780 CONTINUE 

DO 790 IMODE=l, NMODE (IROT) 

^ C0=0 . 0 

w DO 800 ISTAT=1, NSTAT (IROT) 

800 CO=CO+FV( ISTAT, IMODE) **2*WJ (ISTAT) 

^ DO 790 ISTAT=1, NSTAT (IROT) 

W FFV ( IMODE , ISTAT , IROT) =FV ( ISTAT , IMODE) /SQRT (CO) 

SHAPE ( IROT , IMODE , ISTAT , 5 -NTYPE ) =FFV ( IMODE , ISTAT , IROT ) 

IF (NTYPE. NE.l) ZFV ( IMODE , ISTAT , IROT) =FFV (IMODE , ISTAT , IROT) 

«. 79 0 CONTINUE 

— DO 810 ISTAT=1, NSTAT (IROT) 

DO 810 IMODE=l, NMODE (IROT) 

Si 810 WF(ISTAT, IMODE) =WJ (ISTAT) *FFV(IMODE, ISTAT, IROT) 

- *******************************************************************C 

“ DO 820 IMODE=l, NMODE (IROT) 

DO 820 JMODE=l, NMODE (IROT) 
m AT ( IMODE , JMODE ) =0 . 0 

DO 820 ISTAT=1, NSTAT (iROT) 

820 AT (IMODE, JMODE ) =AT ( IMODE, JMODE) + 

„ & WF (ISTAT, JMODE ) *FFV( IMODE, ISTAT, IROT) 

p ™;*** ************************************ ****************************C 
%********* OUTPUT OF TORSIONAL PART ********************************C 
C************* OUTPUT FOR TORSIONAL VIBRATION **********************C 
IF ( NTYPE .EQ. 1 ) THEN 



_ WRITE ( NUNIT (NCOUNT-1) ,*) 7 TORSIONAL MODE SHAPES AND 

& FREQUENCIES 7 
WRITE ( NUNIT (NCOUNT-1) ,*) 

WRITE ( NUNIT (NCOUNT-1) , 1101 ) IROT 
Lo*************** TORSIONAL FREQUENCY ****************************** c 
DO 825 IMODE=l , NMODE ( IROT) 

FREQ ( IROT , IMODE , 5 -NTYPE ) =FR ( IMODE , IROT ) 

825 CONTINUE 

WRITE ( NUNIT (NCOUNT-1) , * ) ' MODE SHAPES ' 

WRITE ( NUNIT (NCOUNT-1) , * ) ( IMODE, IMODE=l, NMODE (IROT) ) 

WRITE ( NUNIT (NCOUNT-1) ,* ) ' TORSIONAL FREQUENCY ( HZ ) ' 

WRITE ( NUNIT (NCOUNT-1) , * ) ( FREQ (IROT, IMODE, 4) /60 . 0 , 

& IMODE=l, NMODE (IROT) ) 

DO 830 ISTAT=1 , NSTAT (IROT) 

— WRITE ( NUNIT (NCOUNT-1) ,* ) ISTAT, (SHAPE (IROT, IMODE , ISTAT, 4 ) , 

& IMODE=l, NMODE (IROT) ) 

830 CONTINUE 
ENDIF 

-O************ OUTPUT FOR Z -DIRECTION VIBRATION ******************************0 
IF ( NTYPE .EQ. 2 ) THEN 

WRITE ( NUNIT (NCOUNT-1) , *) ' Z -DIRECTION MODE SHAPES 7 
WRITE ( NUNIT (NCOUNT-1) ,2101 ) IROT 

~ 0* * ************* Z-DIRECTION FREQUENCY ******************************0 
DO 835 IMODE=l, NMODE (IROT) 

FREQ ( IROT , IMODE , 5 -NTYPE) =ZF ( IMODE , IROT) 

835 CONTINUE 1 

WRITE ( NUNIT (NCOUNT-1) , * ) 7 MODE SHAPES 7 

WRITE ( NUNIT (NCOUNT-1) , * ) ( IMODE, IMODE=l, NMODE (IROT) ) 

W WRITE ( NUNIT (NCOUNT-1) , * ) 7 Z-DIRECTION FREQUENCY ( HZ ) 7 

WRITE ( NUNIT (NCOUNT-1) ,* ) 7 7 , ( FREQ (IROT, IMODE , 3 ) /60 . 0 , 

& IMODE=l, NMODE (IROT) ) 

~ WRITE ( NUNIT (NCOUNT-1) , * ) 

** DO 840 ISTAT=1, NSTAT (IROT) 

WRITE ( NUNIT (NCOUNT-1) , * ) ISTAT, (SHAPE (IROT, IMODE, ISTAT, 3 ) , 

& IMODE=l , NMODE (IROT) ) 

r 840 CONTINUE 
ENDIF 

•'.*********************************************************************C 

61 FORMAT (12, 8F8.3 ) 

^1101 FORMAT ( 7 TORSIONAL FREQUENCY AND ORTHONOMAL MODE ** IROT= 7 ,I2) 

1102 FORMAT (3X, 7 IMODE= 7 , 13 t 30X, 7 FREQUENCY= 7 , 2F10 . 3 ) 

=-1103 FORMAT (2X, 7 ISTAT= 7 , 14X, 7 FV(J,I)= 7 , 12X, 7 FFV(J,I) 7 ,8X, 7 TV= 7 ) 

1104 FORMAT (2X,I3,8X,G14.6;8X,F12.6, 6X , G12 . 5) 

— 1105 FORMAT (3 (/2X, 3F16 . 8) ) 

F1109 FORMAT (2X, 15 , 2F13 . 5) 

~2101 FORMAT ( 7 Z-DIRECTION FREQUENCY AND ORTHONOMAL MODE * IROT= 7 ,I2) 

2102 FORMAT (3X, 7 IMODE= 7 ,I3,30X, 7 FREQUENCY= 7 , F10 . 3 ) 

72103 FORMAT (2X, 7 ISTAT= 7 , 14X, 7 ZV ( J , I ) = 7 , 12X, 7 ZFV ( J , I) 7 , 8X, 7 FZ= 7 ) 

—2109 FORMAT (2X, 15 , F13 . 5 , 2E15 . 8) 

"** CARD FOR THE END OF THE TORSIONAL AND Z-DIRECTION PART **C 
NTYPE= NTYPE + 1 

— IF ( NTYPE-1 .EQ. 1 ) GOTO 700 

1 ******************* CARD FOR THE END OF THE SUBROUTINE ***********C 

*C* CHECK IF THERE ARE MORE ROTORS TO CALCULATE THE MODE SHAPES ****C 
IROT=IROT+ 1 i 

IF ( IROT .LE. NROT ) GOTO 10 

c* ***************************************************************** c 
-******** INITIALIZATION NC9N_BEAR TABLE *****************************0 



DO 201 IROT=l , NROT 

DO 201 I BEAR= 1 , NBEAR ( IROT ) 

DO 201 JROT= 1 , NROT 
DO 201 JBEAR=1, NBEAR (JROT) 
NCON_BEAR ( IROT , I BEAR, JROT , JBEAR) =0 
201 CONTINUE 


;******************************************c 
C* INPUT FOR BEAR CONNECTIONS *************C 


C* 


C* 

LE 

'* 

Ti i* 
c* 

i* 

.* 

't* 

C* 

:* 

wJ* 

c* 

c* 

0 

-j* 


*****************************************c 

READ (NUNIT (1) , *) 

READ (NUNIT (1) , *) JCHECK4 

JCHECK4 IS AN INTEGER SHOWING THE NUMBER OF LINES IN A BEARING CONNECTION TAB 

THE BEARING CONNECTION TABLE IS IN THE FORM: 

IROT IBEAR JROT JBEAR ITYPE_CON 
IROT - I-TH ROTOR INDEX 

IBEAR - I-TH BEARING OF I-TH ROTOR INDEX 

JROT - J-TH ROTOR INDEX 

JBEAR - J-TH BEARING OF J-TH ROTOR INDEX 

ITYPE_CON - TYPE OF CONNECTION BETWEEN (IROT, IBEAR) AND (JROTOR, JBEAR) 

IF ITYPE_CON =0 THEN BEARING IS CONNECTED TO THE GROUND, AND 
THERE IS NO CONNECTION TO ANOTHER BEARING OR TO CASING. 

IF ITYPE_CON ! =0 IT SHOWS THE TYPE OF BEARING STIFFNESS USED 

IN THE CALCULATION OF THE BEARING FORCE BETWEEN (IROT, IBEAR) AND (JROTOR, JBEA 
THE TYPE OF BEARING STIFFNESS IS THE NUMBER OF THE BEARING STIFFNESS TABLE. 


WRITE (NUNIT (NCOUNT) , *) 

WRITE (NUNIT (NCOUNT) ,*) ' BEARINGS CONNECTION TABLE ' 

WRITE (NUNIT (NCOUNT) ,*) ' IROT IBEAR JROT JBEAR 

& ITYPE_CON ' 

IF ( JCHECK4 .EQ. 0 ) THEN 

WRITE (NUNIT (NCOUNT) ,*) ' ALL BEARINGS ARE CONNECTED TO 
& GROUND ’ 

ENDIF 1 


DO 203 JJ=1 , JCHECK4 

IROT, IBEAR CONNECTED TO JROT, JBEAR, TYPE OF CONNECTION 
p READ (NUNIT (1) , *) IROT, IBEAR, JROT, JBEAR, 

w & ITYPE_CON 

NCON_BEAR ( IROT , IBEAR , JROT , JBEAR) =ITYPE_CON 
NCON_BEAR( JROT, JBEAR, IROT, IBEAR) =ITYPE_CON 

~ WRITE (NUNIT (NCOUNT) , *) IROT, ' ' , IBEAR, ' ',JROT,' ' , JBEAR, ' ', 

& ITYPE_CON 

F :1 WRITE ( NUNIT ( NCOUNT ) , * ) 

L203 CONTINUE 

IF ( NSWITCH ( 1) .EQ. 0 ) GOTO 1000 

r t 

•x^** **************************************** ********* t ************* ****c 
c************************* END OF THE MAIN BODY OF LAT ***************0 

^_* **************************************************************** ****Q 

E* WE JUST FOUND ABOVE THE PLANAR UNDAMPED NATURAL FREQUENCIES AND *C 

C* MODE SHAPES . NOW WE PROCEED FOR THE CALCULATIONS OF THE DAMPED *C 

c* MODE SHAPES AND FREQUENCIES WHICH ARE NEEDED FOR THE STABILITY *C 

F* ANALYSIS OF THE PROBLEM. HERE WE JUST CALCULATE A GLOBAL MATRIX *C 

FOR THE EIGENVALUE PROBLEM . USING THIS GLOBAL MATRIX YOU CAN FIND *C 
C* THE EIGENVALUES AND MODE SHAPES BY ANY OTHER STANDART PROGRAM. *C 

.T ; * THIS PART IS ADDED TO CREATE A GLOBAL MATRIX FOR THE EIGENVALUE *C 

SOLVER PROGRAM. *C 

********************************************************************Q 

DO 950 IROT=l , NROT 
H NSPEED= NMODE ( IROT) 



DO 911 1=1 , NSPEED 



DO 911 J=1 , NSPEED 
KBAR (I,«J)=0.0 
CBAR (I , J) =0 . 0 
911 CONTINUE 

DO 900 IBEAR= 1 , NBEAR ( IROT ) 

ISTAT= LBEAR ( IROT , I BEAR ) 

DO 901 ISPEED=1, NSPEED 

FMODE ( IROT , IS PEED , 1 ) = SHAPE ( IROT , I SPEED , ISTAT , 1) 
FMODE ( IROT , ISPEED , 2 ) = SHAPE ( IROT , ISPEED , ISTAT , 2 ) 
FMODE (IROT, ISPEED, 3 ) = FMODE (IROT, ISPEED, 1) 

FMODE (IROT, ISPEED, 4) = FMODE (IROT, ISPEED, 2) 

901 CONTINUE 


— DO 910 ISPEED=1, NSPEED 
DO 910 JSPEED=1, NSPEED 
SBARK= 0.0 

SBARC= 0 . 0 
DO 902 1=1,2 
DO 902 J=1 , 2 
SBARK= SBARK + 

_ & FMODE (IROT, ISPEED, I ) *STIF_BEAR( IROT, IBEAR, I, J) 

& * FMODE (IROT, JSPEED, J) 

S BARC= S BARC + 

& FMODE (IROT, ISPEED, I ) *DAMP_BEAR( IROT, IBEAR, I, J) 

— & *FMODE(IROT, JSPEED, J) 

902 CONTINUE 

KBAR (ISPEED, JSPEED) = KBAR (ISPEED, JSPEED) + SBARK 
** CBAR (ISPEED, JSPEED) = CBAR (ISPEED, JSPEED) + SBARC 

*910 CONTINUE 
900 CONTINUE 



WRITE (56 , * ) ( (CBAR (I, J) ,J=1, NSPEED) ,1=1, NSPEED) 


DO 920 1=1, 2*NSPEED 
r- DO 920 J=l, 2*NSPEED 

~ AG LOB ( IROT , I , J ) = 0.0 

920 CONTINUE 

:** FORM IDENTITY MATRIX BLOCK *C 
DO 925 1=1, NSPEED 
w DO 925 J=NSPEED+1,2*NSPEED 

IF ( I .EQ. J-NSPEED ) AGLOB ( IROT , I , J ) = 1.0 
=925 CONTINUE 

_ DO 930 I=NSPEED+1, 2*NSPEED 

DO 930 J=l, NSPEED 

AGLOB ( IROT , I , J ) = -KBAR (I -NS PEED, J) 

~ IF( J .EQ. I-NSPEED ) 

~ & AGLOB ( IROT , I , J ) = AGLOB ( IROT, I , J) - 

& (FREQ ( IROT, J, l)/9.554)**2.0 
fe930 CONTINUE 

g; DO 935 I=NSPEED+1 , 2*NSPEED 

DO 93 5 J=NSPEED+1 , 2 *NSPEED 
AGLOB ( IROT , I , J) = -CBAR (I-NSPEED , J-NSPEED) 

^ 935 CONTINUE 

F ' 

WRITE ( 30+IROT , *) 7 MATRIX 2*NSPEEDx2*NSPEED ' 

_ WRITE (30+IROT, *) 2*NSPEED, 2*NSPEED 

| WRITE (30+IROT, *) ' GLOBAL MATRIX FOR IROT 7 , IROT 


WRITE ( 3 O+IROT ,940) 

*4 & ( (AGLOB (IROT, I, J) ,J=1, 2 *NSPEED) ,1=1,2 *NSPEED) 

-340 FORMAT ( 6 (Gil. 3, IX) ) 


950 CONTINUE 



x:* ************************************************************* *C 
C******* START CALCULATIONS SUPER-GLOBAL MATRIX **************0 
. ************************************************************** *c 
DO 921 1=1 , 2*NSPEED*NR0T 
” DO 921 J=1 , 2 *NSPEED*NROT 

ASUPER ( I , J ) = 0.0 
921 CONTINUE 

~ DO 952 IROT=l , NROT 

DO 952 1=1 , 2*NSPEED 
DO 952 J=1 , 2*NSPEED 
— 11= (IROT-1) *2*NSPEED + I 

JJ= (IROT-1) *2*NSPEED + J 
ASUPER (II, JJ)= AGLOB ( IROT , I , J) 

952 CONTINUE 


DO 955 IROT=l , NROT 
! ’ DO 956 JROT=l , NROT 

w IF ( IROT .EQ. JROT ) GOTO 956 

;********************************C 

IF ( NMODE (IROT) - NMODE (JROT) .NE. 0 ) THEN 

WRITE (NUNIT (NCOUNT) , *) 'WARNING: ' , IROT , 'AND ',JROT, 

& ' HAVE DIFFERENT NUMBER OF MODES ' 

• WRITE (NUNIT (NCOUNT) , *) 'A-SUPER WILL NOT BE CALCULATED ' 

~ WRITE (NUNIT (NCOUNT) , *) 'CHANGE NUMBER OF MODE SHAPES ' 

STOP 

ENDIF 

■ ************************* *******c 

DO 958 1=1 , NSPEED 
DO 958 J=l, NSPEED 
S K2BAR(I , J) =0 . 0 

C2BAR( I , J) =0 . 0 
«_958 CONTINUE 

S DO 957 IBEAR= 1 , NBEAR ( IROT) 

DO 959 JBEAR= 1 , NBEAR (JROT) 

IF( NCON_BEAR ( IROT, IBEAR, JROT, JBEAR) .NE. 1 ) GOTO 959 
r 5 NSPEED= NMODE (IROT) 

^ ISTAT= LBEAR( IROT, IBEAR) 

JSTAT= LBEAR( JROT, JBEAR) 
n- DO 960 ISPEED=1, NSPEED 

^ FMODE(IROT, ISPEED, 1) = SHAPE (IROT, ISPEED, ISTAT, 1) 

FMODE ( IROT , ISPEED , 2 ) = SHAPE ( IROT , ISPEED , ISTAT , 2 ) 

^ FMODE ( IROT , ISPEED , 3 ) = FMODE ( IROT , ISPEED , 1) 

§§ FMODE (IROT, ISPEED, 4) = FMODE ( IROT , ISPEED, 2 ) 

E FMODE (JROT, ISPEED, 1)= SHAPE (JROT, ISPEED, JSTAT, 1) 

FMODE (JROT, ISPEED, 2) = SHAPE (JROT , ISPEED, JSTAT, 2 ) 
m FMODE (JROT, ISPEED, 3) = FMODE (JROT , ISPEED, 1) 

^ FMODE (JROT, ISPEED, 4) = FMODE (JROT , ISPEED, 2 ) 

'-'960 CONTINUE 

^ DO 961 ISPEED=1, NSPEED 

DO 961 JSPEED=1, NSPEED 

DO 962 1=1,2 

^ SBARK= 0.0 

iri SBARC= 0.0 

— DO 963 J=l, 2 

SBARK= SBARK + 

& ( STIF_BEAR ( IROT , IBEAR, I , J) +STIF_AV ( IROT , IBEAR, I , J) ) 
fe ; & * FMODE (JROT, JSPEED, J) 

SBARC= SBARC + DAMP_BEAR ( IROT , IBEAR, I , J) •* FMODE (JROT , JSPEED, J) 
963 CONTINUE 

f K2BAR (ISPEED, JSPEED) = K2BAR(ISPEED, JSPEED) + 



_ & FMODE (IROT, ISPEED, I) *SBARK 

C2BAR ( IS PEED , JS PEED) = C2BAR ( ISPEED, JSPEED) + 
& FMODE (IROT, ISPEED, I) *SBARC 
962 CONTINUE 
-961 CONTINUE 


959 CONTINUE 
~ 957 CONTINUE 

:****c 

_ DO 964 I=IROT*2*NSPEED, IROT*2*NSPEED+NSPEED-l 

DO 964 J=JROT*2*NSPEED, JROT*2*NSPEED+NSPEED-l 
ASUPER(I-NSPEED+1 , J-NSPEED+1) = 

& C2BAR ( I-IROT*2*NSPEED+l , J-JROT*2*NSPEED+l) 

^ ASUPER ( I-NSPEED+1 , J-2 *NSPEED+1 ) = 

& K2 BAR ( I - IROT * 2 *NSPEED+ 1 , J- JROT* 2 *NSPEED+ 1 ) 


j 964 CONTINUE 
956 CONTINUE 

955 CONTINUE ' 

w WRITE (40, *) 7 MATRIX 2 *NROT*NSPEEDx2 *NROT*NSPEED 7 

WRITE (40,*) 2*NR0T*NSPEED, 2 *NROT*NSPEED 
WRITE (40,*) 7 GLOBAL MATRIX FOR NROT 7 ,NROT 

“ WRITE (40 , 967 ) 

_ & ( (ASUPER(I, J) , J=l,2*NROT*NSPEED) , 1=1 , 2*NROT*NSPEED) 

; 967 FORMAT ( 6(G10.3,1X) ) 


-1000 CONTINUE 
Z: RETURN 

— END 

H SUBROUTINE SBC (NSTATX, NBC, OM, A, GJ , WJ , TK) 

j ***************************** ******************±***** 1 '*******'}:****** 

”c * THIS IS TO SET BOUNDARY CONDITIONS FOR MATRIX TRANSFORM METHOD * 
r ******************************************************************** 

§1 INTEGER NBC 

w REAL D(2,2) ,B(2,2) ,C(2,2) 

REAL GJ (NSTATX) ,WJ( NSTATX ) , TK (NSTATX) 

C (2 , 2 ) =1 . 0 
B ( 1 , 1) =1 . 0 
w B ( 1 , 2 ) =0 . 0 

B (2 , 1) =0 . 0 
ri B (2 ,2) =1.0 

DO 20 1=1, NSTATX 

C(l,l)=1.0+ (TK (I) -WJ (I) *OM**2 ) *GJ (I) 

C ( 1 , 2 ) =TK ( I ) -WJ ( I ) *OM**2 
5-i C (2 , 1) =1 . 0*GJ (I) 

— D(1,1)=C(1,1)*B(1,1)+C(1,2)*B(2,1) 

D(1 , 2 ) =C (1,1)*B(1,2) +C (1,2)*B(2,2) 

ej D ( 2 , 1) =C (2,1)*B(1,1) +C (2,2)*B(2,1) 

D (2 , 2 ) =C (2,1)*B(1,2) +C (2,2)*B(2,2) 

" 6 * B ( 1 / 1) =D( 1 , 1) 

B ( 1 , 2 ) =D (1,2) 
e , B(2 , 1) =D(2 , 1) 

^ B ( 2 , 2 ) =D (2,2) 

"~20 CONTINUE 

NBC=2 

s' FREE FIXED -> NBC=1 

— IF ( NBC .EQ. 1 ) A=B (2,2) 

C FREE FREE -> NBC=2 

= IF ( NBC .EQ. 2 ) A=B (1,2) 



FIXED FIXED -> NBC=3 

IF ( NBC .EQ. 3 ) A=B (2,1) 

FIXED FREE -> NBC=4 

IF ( NBC .EQ. 4 ) A=B(1,1) 

RETURN 

END 

****************************************************************** *c 


